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ABSTRACT 


This  study  investigated  a  (IR)  face  recognition  system  using  an  uncooled  IR  cam¬ 
era.  A  computer-based  image  collection  set-up  was  designed  and  used  to  create  a  small 
database  of  420  facial  images,  from  14  volunteers.  Manual  and  automated  facial  image 
cropping  routines  were  implemented.  Two  linear  approaches  (PCA  and  LDA)  for  the 
dataset  dimension  reduction  and  classification  were  implemented  and  their  resulting  clas¬ 
sification  perfonnances  compared.  Results  show  that  the  best  PCA -based  average  classi¬ 
fication  performance  is  equal  to  92.22%  while  the  LDA-based  classification  performance 
is  equal  to  99.40%.  These  results  successfully  show  that  an  uncooled  IR  camera  may  be 
used  to  discriminate  between  individual  subjects  obtained  from  a  small  database  collected 
under  a  very  controlled  environment. 
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EXECUTIVE  SUMMARY 


This  study  investigated  a  (IR)  face  recognition  system  using  an  uncooled  IR  cam¬ 
era.  A  computer-based  image  collection  set-up  was  designed  and  used  to  collect  a  data¬ 
base.  The  database  has  420  facial  images,  from  14  volunteers  (13  adult  males  and  1  adult 
female)  with  30  pictures  per  person.  All  subjects  are  without  glasses.  Three  different  fa¬ 
cial  expression  sets  were  collected:  first,  the  individuals  wore  a  neutral  expression.  The 
second  set  was  constituted  of  the  same  subjects  with  a  “smiling”  expression.  The  third  set 
was  collected  with  same  individuals  pronouncing  the  vowel  “u”.  Study  variables  were 
reduced  by  using  a  controlled  environment  with  the  same  background  for  each  picture, 
with  the  person-camera  distance  fixed,  and  by  restricting  the  pictures  to  frontal  facial  im¬ 
ages  while  allowing  a  vertical  and  horizontal  angle  freedom  of  10°. 

The  uncooled  infrared  camera  used  for  imaging  has  160  x  120  pixels.  Manual  and 
automated  facial  image  cropping  routines  were  implemented  with  a  fixed  size  of  60  x  45 
pixels.  Two  linear  approaches,  (a)  Principal  Component  Analysis  (PCA)  and  (b)  Linear 
Discriminant  Analysis  (LDA)  for  the  dataset  dimension  reduction  and  classification  were 
implemented  and  their  resulting  classification  perfonnances  compared.  A  minimum  dis¬ 
tance  classifier  was  selected  to  evaluate  the  classification  perfonnances.  The  overall  sys¬ 
tem  performance  was  evaluated  with  a  cross-validation  scheme  using  60%  of  the  pictures 
for  training  and  40%  of  the  pictures  for  testing,  with  1,000  repetitions. 

Results  show  that  the  best  PCA-based  overall  classification  perfonnance 
(92.22%)  is  obtained  when  selecting  the  top  40  eigenvectors,  while  excluding  the  first 
three  top  eigenvectors.  The  LDA-based  approach  performed  better,  with  an  overall  classi¬ 
fication  performance  equal  to  99.40%,  as  expected  from  the  scheme  definition.  Results 
obtained  in  this  study  successfully  show  that  an  uncooled  IR  camera  may  be  useful  to 
discriminate  between  individual  subjects  obtained  from  a  small  database  collected  under 
a  very  controlled  environment. 
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I.  INTRODUCTION 


Infrared  (IR)  imaging  offers  advantages  over  visible  imagery,  the  main  one  is  its 
invariance  to  changes  due  to  illumination  because  the  generated  images  do  not  depend  on 
reflected  light  but  on  radiations  produced  by  the  body.  The  best  IR  image  resolution  is 
produced  using  cooled  camera  equipment,  where  liquid  nitrogen  is  used  to  reduce  the 
temperature  of  the  sensors.  IR  face  recognition  using  cooled  IR  imagery  has  been  used 
for  several  years;  however,  the  high  associated  equipment  cost  has  slowed  down  its  use 
significantly.  However,  recent  technological  developments  have  resulted  in  uncooled  IR 
camera  set-ups  with  resolution  approaching  that  of  uncooled  cameras  at  a  fraction  of  the 
cost.  This  study  investigates  the  design  of  a  (IR)  face  recognition  system  using  an  un¬ 
cooled  IR  camera. 

First,  Chapter  II  describes  the  applications  of  biometrics  to  everyday  life.  Chapter 
III  presents  infrared  imaging  concepts  and  the  hardware  set-up  selected  in  our  study  to 
conduct  the  experiments.  We  also  describe  in  this  Chapter  the  image  grabbing  procedure 
and  the  nomenclature  selected  for  the  database  collection.  Chapter  IV  presents  the  overall 
face  recognition  algorithms  implemented  for  this  stud.  First,  we  describe  two  face- 
extracting  algorithms,  which  allow  isolating  the  image  face-only  portions.  Next,  we  pre¬ 
sent  the  two  dimension  reduction  schemes  considered  in  this  work,  Principal  Component 
Analysis  (PCA)-based  and  the  Linear  Discriminant  Analysis  (LDA).  The  minimum  dis¬ 
tance  classifier  approach  followed  to  make  class  decisions.  Chapter  V  presents  the  recog¬ 
nition  performances  obtained  on  our  database  and  discusses  the  cross-validation  approach 
implemented  to  estimate  the  performance.  Conclusions  and  recommendations  for  future 
study  are  presented  in  Chapter  VI.  Finally,  all  code  implemented  during  the  study  is  in¬ 
cluded  in  the  Appendix. 
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II.  BIOMETRICS 


This  section  presents  the  basic  human  biometrics  used  today,  each  with  their  ad¬ 
vantages  and  drawbacks. 


A.  INTRODUCTION 

The  September  11th  terrorist  attack  clearly  demonstrated  that  security  in  a  broad 
spectrum  of  areas  needs  to  be  intensified  and  that  the  positive  identification  of  individuals 
is  highly  desirable  in  many  instances.  For  example,  identification  can  prevent  the  unau¬ 
thorized  access  to  confidential  information  or  physical  facilities  and  can  help  prevent  ter¬ 
rorism  and  crime. 

Biometrics,  can  certainly  increase  the  probabilities  of  positive  identification  by 
using  physical  and  behavioral  characteristics  to  identify  persons  or  verify  identities.  Some 
of  the  physical  and  behavioral  characteristics  used  in  biometrics  are  the  following: 

•  Physical  Characteristics:  chemical  composition  of  body  odor;  facial  features 
and  infrared  thermal  emissions;  features  of  the  retina  and  iris;  fingerprints; 
hand  geometry;  wrist  and  hand  veins. 

•  Behavioral  Characteristics:  handwritten  signature;  keystrokes  dynamics  and 
voiceprint. 

Despite  widespread  media  coverage  of  biometrics  since  the  September  11th  terrorist  at¬ 
tacks,  a  national  survey  conducted  by  SEARCH,  the  National  Consortium  for  Justice  In¬ 
formation  and  Statistics,  showed  that  only  half  of  the  general  public  is  aware  of  such 
technologies.  A  public  opinion  poll  was  conducted  between  September  18th  and  30th, 
2001,  shortly  after  the  terrorist  attack  and  again  between  August  15th  and  18th,  2002.  Ac¬ 
cording  to  these  surveys,  public  support  for  using  biometrics  for  anti-terrorism  measures 
or  crime  prevention  declined  from  86%  in  2001  to  80%  in  2002  [1].  In  addition,  it  is  in¬ 
teresting  to  note  that  these  facts  indicate  that  biometric  identification  systems  remain  rela¬ 
tively  unknown,  yet  are  generally  supported.  However,  the  U.S.  Congress  has  recently 

passed  several  laws  that  are  likely  to  make  biometrics  far  more  common.  These  laws  are 
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only  likely  to  be  the  inception  of  an  inevitable  trend  to  employ  Biometrics  for  the  pub¬ 
lic’s  safety.  For  example: 


•  The  USA  Patriot  Act,  passed  in  November  2001,  requires  the  Federal  Gov¬ 
ernment  to  develop  technology  standards  for  verifying  visa  applicants  identiti- 
ties  to  ensure  that  a  specific  individual  has  not  received  a  visa  under  a  differ¬ 
ent  name; 

•  The  Enhanced  Border  Security  and  Visa  Entry  Reform  Act,  which  President 
Bush  signed  in  May  2002,  requires  that  visas  and  other  travel  documents  in¬ 
clude  biometric  identifiers  by  October  26,  2004; 

•  The  Aviation  and  Transportation  Security  Act,  signed  in  November  2001,  au¬ 
thorizes  the  use  of  biometrics  or  other  technologies  to  verify  identities  of  those 
entering  secure  airport  areas  [2], 


Furthermore,  additional  funding  is  being  allocated  to  develop  technologies  and 
policies  focused  on  combating  cyberterrorism.  For  example,  the  Defense  Department 
awarded  Carnegie  Mellon  University  a  five-year,  $35.5  million  grant  in  that  area  in  Sep¬ 
tember  2002.  According  to  Praddeep  Khosla,  director  of  the  Center  for  Computer  and 
Communications  Security,  the  research  will  focus  on  four  key  areas,  one  of  which  being 
computer  secured  access,  and  biometrics-based  secured  access  to  various  devices,  and 
public  and  private  facilities  [3]. 

The  International  Biometric  Association  estimates  that  sales  in  the  biometric  in¬ 
dustry  reached  $170  million  in  2001,  a  70%  increase  over  the  previous  year.  Furthermore, 
it  also  predicts  that  sales  will  rise  to  $1  billion  in  2004  [4].  According  to  findings  from 
the  international  market  consulting  company,  Frost  &  Sullivan,  the  total  biometric  market 
is  projected  to  reach  $2.05  billion  in  2006  [5]. 

Biometrics  usually  attempt  to  solve  two  problems,  identification  and  verification. 
Identification  is  the  capability  to  properly  identify  a  person  from  all  those  whose  biomet¬ 
ric  measurements  have  been  previously  collected  and  stored  in  a  database.  Verification  is 
authenticating  a  person’s  identity  by  comparing  the  biometrics  stored  on  a  database. 
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Some  biometric  properties  are  desirable  but  are  not  completely  attainable  with  to¬ 
day’s  technology.  There  are  biometrics  that  attain  completely  some  but  not  all  properties. 
Note  that  some  very  accurate  biometrics  may  not  be  acceptable,  which  reduces  their 
range  of  applications.  Conversely,  some  less  accurate  biometrics  may  be  used  as  they  are 
easier  to  implement  and  are  more  widely  accepted.  For  example,  a  gym  academy  could 
use  a  biometric  system  to  restrict  its  access.  It  is  possible  to  use  a  lower  performing  bio¬ 
metrics  that  is  more  acceptable  for  such  an  application,  where  granting  incorrect  access 
has  limited  consequences.  However,  such  a  system  would  not  likely  be  employed  in 
places  that  require  a  more  controlled  access  such  as  nuclear  power  plants. 


B.  BIOMETRICS  PROPERTIES 

Any  human  physiological  or  behavioral  characteristic  can  be  used  in  biometrics, 
provided  it  has  the  following  desirable  properties: 

•  Universality:  all  individuals  requiring  identification  or  verification  must  have 
their  specific  biometric  characteristic  measured; 

•  Uniqueness:  no  two  persons  can  possess  identical  biometric  characteristics; 

•  Pennanence:  the  characteristic  must  be  invariant  with  time; 

•  Collectability:  the  characteristic  can  be  measured  quantitatively. 

Other  requirements  are  also  important,  such  as: 

•  Performance:  refers  to  how  accurately  the  biometric  data  can  be  identified  or 
verified; 

•  Acceptability:  indicates  to  what  extent  people  are  willing  to  accept  the  intru¬ 
siveness  of  the  biometric  system  without  feeling  their  privacy  or  personal 
freedoms  are  being  breached; 

•  Circumvention:  refers  to  how  easily  the  biometrics  system  can  be  deceived. 
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C.  BIOMETRIC  TYPES 


Biometrics  most  commonly  studied  are  the  face,  fingerprint,  hand  geometry,  iris, 
retinal  pattern,  signature,  voiceprint  and  thermograms.  Note,  however,  that  no  single 
biometric  is  universal.  For  example,  some  individual  may  have  a  finger  skin  too  dry  to 
gather  fingerprints,  as  may  be  the  case  for  a  hairdresser  routinely  working  with  harsh 
chemicals  [4].  Other  individuals  may  be  mute  or  amputees.  The  most  common  biomet¬ 
rics  are  described  as  follows: 

1.  Voice 

The  voice  is  a  unique  human  characteristic.  However,  it  is  usually  degraded  in 
quality  by  microphone,  communication  channel  or  analog  digital  converters  used  to  col¬ 
lect  and  to  transmit  information.  However,  voice  acquisition  schemes  are  non-intrusive 
and  many  companies,  including  Bell  Canada  and  Visa  International,  have  tested  voice- 
print  technology.  Visa  could  eventually  use  this  technology  to  identify  credit  card  holders 
by  having  them  speak  into  their  PCs  microphones  to  verify  online  purchases. 

Visa  has  been  testing  the  software  Vocent’s  Voice  Secure  from  Vocent  Solutions 
Inc.  of  Mountain  View,  California.  Its  initial  usage  has  been  used  solely  within  the  firm, 
but  the  company  is  also  considering  the  idea  of  deploying  the  software  as  a  customer¬ 
facing  application  next  year  [6]. 

Bell  Canada  selected  Nuance  Verifier  3.0  voice  authentication  software  as  a  secu¬ 
rity  solution  for  the  Speech  Enabled  Field  Access  System.  This  voiceprint  authentication 
technology  allows  5,000  Bell  Canada  technicians  to  access  the  dispatch  system,  and  up¬ 
date  customer  installation  and  repair  orders  securely  from  any  telephone  [7]. 

However,  voiceprints  have  their  own  limitations  as  reproductions  of  a  recorded 
voice  can  circumvent  a  voice  authentication  system  for  a  remote  unattended  application. 
In  such  cases,  additional  security  may  be  obtained  by  varying  the  spoken  sentence. 
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2.  Face 


The  face  is  the  most  acceptable  biometrics  because  it  is  one  of  the  most  common 
methods  of  identification,  which  humans  use  daily.  Furthermore,  the  method  for  acquir¬ 
ing  samples  for  verification  is  non-intrusive.  However,  developing  face  recognition  tech¬ 
niques  that  can  allow  for  the  effects  of  aging,  facial  expressions,  slight  variations  in  imag¬ 
ing  and  photographing  variations  in  the  facial  pose  is  quite  challenging.  A  facial- 
recognition  system  tested  in  2002  at  a  Palm  Beach  airport  in  Florida  illustrates  the  sys¬ 
tem’s  limitations.  The  system  failed  to  match  airport  employs  with  their  digital  photos 
53%  of  the  time  [8].  Nevertheless,  in  March  2002  the  U.S.  Army  announced  that  its  mili¬ 
tary  police  officers  have  successfully  tested  facial-recognition  technology  to  aid  them  in 
their  duties.  The  cornerstone  of  this  mobile  security  system  is  an  eyeglass-mounted  wear¬ 
able  camera  and  display  device. 

According  to  Mark  Chandler,  a  physical  scientist  at  the  Anny’s  Soldier  Systems 
Center  in  Natick,  Massachusetts,  the  Army’s  development  of  this  mobile  facial  recogni¬ 
tion  system  has  been  a  four-year,  $12  million  effort.  The  U.S.  Army  intends  to  integrate  it 
with  other  technologies,  and  the  U.S.  Navy  is  investigating  its  applications  to  control  ac¬ 
cess  to  flight  lines  [9]. 


3.  Fingerprint 

Fingerprints  are  believed  to  be  unique  to  each  person  (and  each  finger).  As  a  re¬ 
sult,  fingerprinting  is  one  of  the  most  mature  biometric  technologies  used  in  forensic  sci¬ 
ence  worldwide.  Fingerprint  images  are  captured,  typically,  in  one  of  two  ways,  by  scan¬ 
ning  an  inked  impression  of  a  finger  or  by  using  a  live-scan  fingerprint  scanner  [10]. 

The  FBI  inaugurated  a  fully  operational  Integrated  Automated  Fingerprint  Identi¬ 
fication  System  (IAFIS)  in  1999  with  development  and  implementation  costs  around 
$640  million.  The  IAFIS  has  its  roots  in  the  1960s  and  1970s  when  the  FBI  began  inves¬ 
tigating  the  feasibility  of  automating  the  fingerprint  identification  process  [11,12].  This 
system  compares  fingerprints  against  those  contained  in  its  vast  database  and  can  respond 
within  two  hours.  In  the  past,  this  process  often  took  weeks  to  complete  using  a  semi- 
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automated  system.  In  addition  to  criminal  checks,  the  FBI  also  processes  many  civil  fin¬ 
gerprint  submissions.  The  FBI  receives  an  average  of  50,000  fingerprint  submissions 
each  business  day,  where  half  on  those  submissions  are  criminal  and  half  are  civil.  Note 
that  civil  submissions  are  those  for  which  law  requires  an  FBI  fingerprint  check,  such  as 
for  teachers,  child  care  providers,  security  guards  and  other  such  positions  of  trust.  The 
IAFIS  now  provides  responses  to  civil  submissions  within  24  hours,  while  such  requests 
response  times  exceeded  100  days  before.  In  addition,  fingerprint  recognition  systems 
have  also  been  employed  for  everyday  applications,  such  as  for  personal  computers,  gro¬ 
cery  stores,  etc.  For  example,  customers  have  the  option  to  pay  for  their  purchases  by 
placing  their  fingers  on  a  payment  terminal  at  the  checkout  register  at  three  Kroger  stores 
in  College  Station,  Texas  and  West  Seattle  Triftway  [13].  The  Children’s  Hospital  Inc.  in 
Columbus,  Ohio,  uses  fingerprint  technology  to  authenticate  doctors  who  order  tests  and 
prescribe  drugs  [14]. 

Semiconductor  makers  have  developed  chips  that  can  read  a  fingerprint  and  verify 
a  specific  user  identity.  For  example,  chips  made  by  STMicroelectronics  (called  the 
Touch  Chip)  and  Veridicom  have  already  been  installed  on  laptops  from  MicronPC  LLC, 
Samsung  Corp.,  NEC  Corp.  and  Acer  Inc.  [15]. 

4.  Thermograms 

A  thermogram  measures  the  temperature  of  the  body’s  surface,  and  as  a  result,  the 
face  thermogram  may  be  used  as  a  biometric.  Prokoski,  a  researcher  from  Mikos  Ltd., 
showed  that  the  temperature  pattern  on  the  facial  skin  is  due  primarily  to  the  pattern  of 
superficial  blood  vessels,  which  are  those  under  the  skin  but  above  bone  and  muscle  [16]. 
Vessels  transport  warm  blood  throughout  the  body  and  heat  the  skin.  As  a  result,  the  skin 
that  is  directly  above  a  blood  vessel  is  on  the  average  0.1°C  warmer  than  the  adjacent 
skin.  The  typical  human  face  displays  an  apparent  temperature  range  of  about  8°C. 

The  complexity  of  the  resulting  3.1  miles  (five  kilometers)  of  bloods  vessel  in  the 
head  and  face  assures  that  each  person’s  vascular  arrangement  is  irreproducible  and 
hence  unique.  Amazingly,  studies  have  shown  that  even  identical  twins  have  different 

thennogram  features  [16].  Furthermore,  studies  have  shown  that  the  skin  temperature  ex- 
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hibits  several  periodic  fluctuations  and  that  the  mechanisms  and  physiological  basis  un¬ 
derlying  these  fluctuations  are  not  yet  well  understood.  Additionally  experiments  indicate 
that  personality  and  emotional  state  are  factors  that  contribute  to  these  fluctuations  [17, 
18]. 


5.  Retinal  Scan 

The  retinal  vasculature  is  rich  in  structure  and  is  supposed  to  be  a  characteristic  of 
each  individual  and  each  eye.  It  is  claimed  to  be  the  most  secure  biometrics  since  it  is  not 
easy  to  alter  or  replicate.  However,  many  people  are  suspicious  of  retinal  scans  owing  to 
the  fact  that  retinal  vasculature  can  reveal  some  medical  conditions  such  as  hypertension. 
This  fact  presently  stands  in  the  way  of  public  acceptance  of  retinal  scan  based- 
biometrics. 

6.  Iris 

The  visual  texture  of  the  human  iris  is  determined  during  the  embryonic  develop¬ 
ment  and  is  posited  to  be  unique  for  each  person  and  each  eye  [10].  Unlike  retinal  scans, 
which  require  close  contacts  with  the  scanner,  iris-based  recognition  has  been  reported  to 
be  successful  at  distances  of  up  to  18  inches  (46  cm)  [19].  However,  iris  scan  systems 
may  be  disrupted  by  contact  lenses  [20]. 

Iridian  Technologies  Inc.,  produces  the  IrisAcess  2200  system,  which  detects  an 
individual  approaching  the  imager  and  then  captures  an  iris  image  once  the  person’s  eye 
is  from  three  to  10  inches  from  the  mirror  on  the  unit.  This  image  is  digitally  processed 
into  a  512-byte  IrisCode  template  [21]. 

The  Pentagon  tested  iris  recognition  technology  at  its  Athletic  Club,  as  club 
members  may  voluntarily  sign  up  to  test  the  system,  which  involves  capturing  data  from 
a  member’s  identification  card  and  iris  features.  In  addition,  Canada  customs  will  begin 
using  iris  scanners  to  speed  air  travelers  through  the  country’s  busiest  airports  in  the  near 
future.  Kiosks  equipped  with  iris-recognition  devices  in  eight  of  Canada’s  busiest  intema- 


9 


tional  airports  will  allow  some  Canadian  travelers  to  cross  customs  checkpoint  in  30  sec¬ 
onds  or  less  by  confirming  their  identities  with  quick  iris  scans  [22], 

7.  Ear 

Studies  have  shown  that  the  shape  of  the  ear  and  the  structure  of  the  pinna  carti¬ 
laginous  tissue  are  distinctive.  In  addition,  there  is  supporting  evidence  indicating  that 
ears  are  probably  unique  to  each  individual.  A.  Iannarelli  first  developed  an  antropomet- 
ric  ear-based  biometric  identification  technique  in  1949.  This  method  is  based  on  12 
measurements  of  the  ear  structure.  However,  this  method  is  not  well  suited  for  “machine 
vision”  implementations  due  to  the  difficulties  of  locating  the  anatomical  point  that 
serves  as  the  origin  of  the  measurement  system  [19].  Although  this  research  is  not  yet 
conclusive,  studies  have  also  indicated  that  identical  twins  have  different  ear  structures 
(even  though  they  may  look  similar).  [19]. 

8.  Gait 

Gait  refers  to  the  peculiar  way  one  walks  and,  while  it  is  not  claimed  to  be  abso¬ 
lutely  unique  to  each  individual,  it  is  sufficiently  characteristic  to  allow  for  identity  au¬ 
thentication.  Naturally,  gait  is  a  behavioral  biometric  and  may  not  stay  invariant  espe¬ 
cially  over  a  long  period  of  time,  due  to  large  fluctuations  of  body  weight,  major  injuries 
involving  the  joints  or  cerebral  disorders  or  inebriety.  Typically,  gait  features  are  derived 
from  analyzing  video-sequence  footage  of  a  person  walking  and  encompass  several  char¬ 
acteristic  movements  of  each  articulate  joint. 

Gait  recognition  research  has  largely  been  motivated  by  Johansson’s  experiment 
and  the  ability  of  humans  to  perceive  motion  from  Moving  Light  Displays  (MLDs).  In 
these  experiments,  human  subjects  were  able  to  recognize  the  type  of  movement  of  a  per¬ 
son  by  observing  the  2D  motion  pattern  generated  by  light  bulbs  attached  to  the  person. 
Similar  experiments,  later  showed  evidence  that  the  identity  of  a  familiar  person,  as  well 
as  the  gender  of  the  person  might  be  recognized  from  MLDs  [23]. 
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9.  Keystroke  Dynamics 


This  behavioral  biometrics  is  also  not  expected  to  be  unique  to  each  individual, 
but  it  may  offer  sufficient  discriminatory  infonnation  to  permit  for  identity  authentica¬ 
tion.  Keystroke  dynamic  features  are  based  on  time  durations  between  the  keystrokes. 
Some  variants  of  identity  authentication  use  features  based  on  inter-key  delays  as  well  as 
dwell  times — namely,  how  long  a  person  holds  down  a  key. 


10.  DNA 

DNA  (Deoxyribonucleic  Acid)  is  mostly  used  in  forensics  for  identification  and  is 
believed  to  be  unique  for  each  person  (except  for  identical  twins  who  have  identical  DNA 
patterns).  Three  issues  limit  the  utility  of  this  biometric  for  other  applications  [10]: 

•  Contamination  and  Sensitivity:  it  is  easy  to  steal  a  DNA  sample  for  a  fraudu¬ 
lent  or  criminal  purpose; 

•  Automatic  Real-time  Identification  Issues:  presently,  the  technology  for  ge¬ 
netic  matching  for  online  identification  is  lacking.  Most  of  the  human  DNA  is 
identical  for  the  entire  human  species  and  only  some  relatively  small  numbers 
of  specific  locations  on  DNA,  called  polymorphic  loci,  exhibit  individual 
variations.  These  variations  are  manifested  either  in  the  number  of  repetitions 
of  a  block  of  base  sequences  or  in  the  minor  non- functional  perturbations  of 
the  base  sequence.  The  processes  involved  in  DNA-based  personal  identifica¬ 
tion  determine  if  two  DNA  samples  originate  from  the  same  individual  based 
on  the  distinctive  signature  at  one  or  more  polymorphic  loci.  The  major  com¬ 
ponents  of  these  processes  now  exist  in  the  form  of  chemical  methods  requir¬ 
ing  an  expert’s  skills. 

•  Privacy  Issues:  information  about  susceptibilities  of  a  person  to  certain  dis¬ 
eases  could  be  gained  from  the  DNA  pattern,  and  there  is  a  significant  concern 
that  the  unintended  abuse  of  genetic  code  infonnation  may  result  in  discrimi¬ 
nation,  such  as  in  hiring  practices. 
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11.  Signature 


The  way  a  person  signs  his  or  her  name  is  known  to  be  a  characteristic  of  that  in¬ 
dividual.  However,  signatures  are  not  a  foolproof  behavioral  biometric  as  they  evolve 
over  a  period  of  time  and  are  influenced  by  physical  and  emotional  conditions.  Although 
experts  can  discriminate  between  genuine  and  forged  signatures,  modeling  the  invariance 
in  the  signatures  and  automatic  signature  recognition  process  are  challenging.  Two  ap¬ 
proaches  to  signature  verification  exist,  static  and  dynamic.  In  the  static  signature  verifi¬ 
cation,  only  geometric  features  of  the  signature  are  used  for  authentication.  In  dynamic 
signature  verification,  the  signature  dynamic  features  such  as  acceleration,  velocity,  and 
trajectory  profiles  are  also  employed. 

A  related,  and  remarkable,  technology  is  authentication  based  on  the  characteris¬ 
tics  of  the  sound  produced  during  a  signature.  These  acoustic  emissions  are  claimed  to  be 
a  characteristic  of  each  individual  [10]. 

12.  Odor 

It  is  commonly  known  that  body  odors  serve  several  functions  including  commu¬ 
nication,  attracting  mates,  assertion  of  territorial  right,  and  protection  from  a  predator. 
Humans  produce  an  identifiable  odor  that  is  characteristic  of  their  unique  chemical  com¬ 
position.  A  component  of  the  odor  emitted  by  human  (or  any  animal)  is  distinctive  to  a 
particular  individual.  Despite  the  rare  use  of  this  technology  in  biometrics,  there  are  many 
commercial  electronic  “noses”  and  software  that  can  recognize  odors  under  development. 
Generally,  these  common  applications  are  in  food  and  beverage  processing,  environ¬ 
mental  monitoring,  medical  diagnostic,  and  fragrance  development  [24,  25]. 

The  artificial  nose  manufactured  by  SMart  Nose  Inc.,  uses  mass  spectrometry  to 
identify  volatile  organic  components  ranging  from  liquid  to  solid  samples.  The  instru¬ 
ment  is  fully  automated  and  operates  in  two  modes:  a  survey  mode  for  the  setup  of  new 
applications,  and  a  multiple  channel  mode  for  the  specific  and/or  sensitive  analyses  in  the 
day-to-day  use.  The  system  uses  the  algorithm  based  on  the  Principal  Component  Analy¬ 
sis  (PC A)  or  Discriminant  Function  Analysis  (DFA)  [24]. 
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13.  Hand  and  Finger  Geometry 


Some  features  related  to  a  human  hand  such  as  finger  length  are  relatively  invari¬ 
ant  and  peculiar  (although  not  unique)  to  each  individual  [10].  Hand  geometry  readers 
have  been  the  dominant  biometric  technology  for  access  control  and  time  and  attendance 
applications  in  2001,  as  reported  by  Frost  &  Sullivan’s  World  Biometric  Report  2002 
[26].  Since  1991,  HandReaders  from  IR  Recognition  Systems,  have  produced  more  than 
100  million  biometric  verifications,  with  more  than  50,000  produced  on  high  volume 
days  at  the  San  Francisco  International  Airport.  HandReaders  span  the  entire  airport,  se¬ 
curing  more  that  180  doors  and  verifying  the  identity  of  more  than  18,000  employees. 

Over  70,000  HandReaders  are  installed  throughout  the  world  in  a  wide  variety  of 
applications.  The  1996  Olympic  games  used  HandReaders  to  protect  access  to  the  Olym¬ 
pic  Village,  where  more  than  65,000  people  were  enrolled  and  over  1  million  transactions 
were  handled  in  28  days  [26]. 


D.  COMPARISON  AMONG  BIOMETRICS 

Biometrics  may  be  compared  in  relation  to  their  desirable  properties.  Table  1  pro¬ 
vides  a  comparison  of  biometrics  technologies. 
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Table  1.  Biometric  Technologies  Comparison.  The  data  was  obtained  from  [10]. 


Biometrics 

Universality 

Uniqueness 

Permanence 

Collectability 

Performance 

Acceptability 

Circumvention 

Face 

High 

Low 

Medium 

High 

Low 

High 

Low 

Fingerprint 

Medium 

High 

High 

Medium 

High 

Medium 

High 

Hand 

Geometry 

Medium 

Medium 

Medium 

High 

Medium 

Medium 

Medium 

Keystrokes 

Low 

Low 

Low 

Medium 

Low 

Medium 

Medium 

Hand  Vein 

Medium 

Medium 

Medium 

Medium 

Medium 

Medium 

High 

Iris 

High 

High 

High 

Medium 

High 

Low 

High 

Retinal  Scan 

High 

High 

Medium 

Low 

High 

Low 

High 

Signature 

Low 

Low 

Low 

High 

Low 

High 

Low 

Voice  Print 

Medium 

Low 

Low 

Medium 

Low 

High 

Low 

Face  Ther¬ 
mograms 

High 

High 

Low 

High 

Medium 

High 

High 

Odor 

High 

High 

High 

Low 

Low 

Medium 

Low 

DNA 

High 

High 

High 

Low 

High 

Low 

Low 

Gait 

Medium 

Low 

Low 

High 

Low 

High 

Medium 

Ear 

Medium 

Medium 

High 

Medium 

Medium 

High 

Medium 

This  Chapter  provided  an  overview  of  the  use  of  biometrics  today  with  their  own 
advantages  and  drawbacks.  Next,  we  discuss  the  main  concepts  behind  infrared  imaging. 
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III.  INFRARED  IMAGING 


This  chapter  presents  the  infrared  imaging  system  used.  First,  we  address  relation¬ 
ships  between  the  temperature  of  a  source  and  the  infrared  radiation.  Second,  we  present 
the  camera  and  physical  set-up  used  to  obtain  the  images. 


A.  INFRARED  SPECTRUM 

The  astronomer  Sir  William  Herschel  (1738  -  1822)  was  the  first  to  detect  infra¬ 
red  (IR)  radiation  in  1800.  IR  radiation  covers  0.7  to  1000  pm  in  the  electromagnetic 
spectrum.  According  with  Hecht,  the  IR  spectrum  is  subdivided  into  four  regions:  the 
near  IR  (0.78  -  3  pm),  the  intermediate  IR  (3  -  6  pm),  the  far  IR  (6  -  15  pm)  and  the  ex¬ 
treme  IR  (15  -  1000  pm)  [27]. 

B.  BLACKBODY 

A  blackbody  is  a  perfect  radiator.  It  radiates  more  power  from  a  surface  area  in  a 
wavelength  interval  than  any  body  can  radiate  at  a  given  temperature.  No  surface  that  is 
in  thermodynamic  equilibrium  can  radiate  more  power,  unless  it  contains  fluorescent  or 
radioactive  materials. 

The  expression  for  the  blackbody  emission  was  first  formally  derived  by  Planck 
in  1900.  The  exitance  for  a  blackbody,  M(A,T),  considering  the  irradiation  in  a  hemi¬ 
sphere  in  front  of  the  source,  is  given  by  [28]: 


M(A,T) 


watt 
cm2  •  pm 


Inhc1 

A\ehc,m  _!]  ’ 


(3.1) 
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where  h  =  6.625xl0'34  Watts-  s2  is  the  Plank’s  constant,  k  =  1.38x10  °3  Watts-s-K  'is 

o 

the  Boltzman’s  constant,  c  =  3x10  m/s  is  the  speed  of  light,  and  T  is  the  temperature  in 
Kelvin.  M(/t,T)  is  usually  expressed  in  Watt  •  cm’  •  pm.  Figure  1  plots  the  exitance  as  a 
function  of  the  wavelength,  for  three  different  temperatures  near  room  temperature 
(300K). 
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Figure  1 .  Blackbody  exitance. 


The  Wien  displacement  law  defines  the  wavelength  that  corresponds  to  the  peak 
of  the  exitance  at  each  temperature.  This  can  be  obtained  by  computing  the  partial  deriva¬ 
tive,  with  respect  to  wavelength,  of  the  exitance  and  setting  it  equal  to  zero,  expressed  as: 

^-0,  (3.2) 

dA 

which  leads  to  /lmax T  =  2898  pm  •  K,  where  Amax  is  the  wavelength  where  the  exitance  peak 

occurs  at  the  temperature  T.  The  maximum  exitance  occurs  at  /L  =  9.66  pm  for  a  black- 

body  source  at  T  =  300  K  (room  temperature),  while  the  maximum  exitance  occurs  at  A. 

=  9.35  pm  for  a  blackbody  source  at  T  =  3 10  K  (human  body  temperature). 
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The  total  blackbody  exitance  is  obtained  by  integrating  the  exitance  expression 
over  a  finite  spectral  region,  which  leads  to: 

^1 

M  =  $M(A,T)dA,  (3.3) 

where  A/  and  A2  are  the  smallest  and  largest  wavelengths  of  the  finite  spectral  region  of 
interest. 

The  total  exitance  from  a  blackbody  at  temperature  T  is  the  integral  over  all  wave¬ 
lengths  values,  and  is  given  by: 

oo 

M(T)  =  \M(A,T)dA,  (3.4) 

0 

which  leads  to  M(T)  =  cieT4,  where  <je  is  the  Stefan-Boltzmann  constant  which  has  an  ap- 

12  2  4 

proximate  value  of  5.67x10'  watt/(cm“-  K  ). 


C.  EMISSIVITY 


The  previous  section  presented  the  equation  of  the  blackbody  exitance.  Recall  that 
the  blackbody  is  a  perfect  radiator  as  it  radiates  the  maximum  power  from  a  surface  area 
in  a  wavelength  interval.  So,  it  also  represents  an  upper  bound  for  the  overall  exitance  of 
a  source. 


The  ratio  between  the  exitance  of  a  real  source  and  the  exitance  of  a  blackbody  is 
defined  as  the  emissivity  s,  which  is  a  function  of  the  wavelength  A  and  the  absolute  tem¬ 
perature  T.  It  is  a  dimensionless  quantity  smaller  than  1  and  defined  as: 


s(A,T)  = 

M(AJ)blackbody 


(3.5) 


Note  that  a  blackbody  source  has  an  emissivity  s  equal  to  1  for  all  wavelengths, 
while  a  graybody  source  has  an  emissivity  less  than  1  for  all  wavelengths. 
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Human  skin  emissions  can  be  approximated  by  a  graybody  source  with  s  equal  to 
0.98  in  the  wavelength  range  between  3  and  14  pm  [29].  For  a  typical  person,  tempera¬ 
ture  across  the  face  varies  about  8°C  [16].  Since  this  variation  is  much  smaller  than  the 
average  absolute  temperature  (Tav)  of  the  face  (310  K),  the  variation  of  the  exitance 
across  the  face  can  be  taken  expressed  as: 

AM  -  4a  AT .  (3.6) 

For  example,  the  face  exitance  changes  by  approximately  1%  when  the  face  tem¬ 
perature  changes  by  1%. 


D.  INFRARED  IMAGES 

1.  IR  Camera 

The  uncooled  IR  camera  selected  for  this  study  is  a  IR-160  model  produced  by  In¬ 
frared  Solutions,  Inc.,  where  the  primary  reason  for  using  an  uncooled  camera  was  to  re¬ 
duce  the  overall  system  cost.  This  uncooled  camera  is  sensitive  to  electromagnetic  radia¬ 
tion  in  the  wavelength  range  from  8  pm  to  14  pm.  It  has  a  160  x  120  pixel  microbolome¬ 
ter  array  to  obtain  the  images  at  a  30  Hz  frame  rate  that  can  be  displayed  on  a  NTSC 
standard  monitor,  or  may  be  transmitted  serially  using  8  bits  per  pixel  via  RS-232.  Mi¬ 
crobolometers  operate  on  the  principle  that  temperature  changes  produced  by  the  ab¬ 
sorbed  radiation  produce  changes  in  the  material  resistance  [28].  As  a  final  result,  an 
electrical  signal  proportional  to  the  change  in  the  electrical  resistance  is  obtained.  The 
bolometer  temperature  coefficient  depends  on  the  sensor  material  and  the  amount  of  re¬ 
sistance  change  produced  by  a  temperature  change.  The  relation  between  the  temperature 
and  the  resistance  may  be  approximated  by  the  linear  equation: 

R  =  R0(l  +  aATd),  (3.7) 

where  a  is  the  temperature  coefficient  of  the  detector  material,  R0  is  the  resistance  at  the 
base  temperature  and  ATci  is  the  change  in  the  detector  temperature.  The  camera  has  a 
germanium  lens  that  allows  for  manual  focusing  of  the  image.  In  addition,  it  is  possible  to 
adjust  the  image  brightness  and  contrast  and  to  select  a  specific  palette  scheme.  Selec- 
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tions  are  saved  in  non-volatile  memory  and  the  camera  initiates  with  the  same  setting  the 
next  time  it  is  turned  on.  The  palette  scheme  used  in  our  study  is  linear,  which  allows  us 
to  adjust  the  gray  color  displayed  automatically.  In  addition,  it  sets  the  output  image  pixel 
that  corresponds  to  the  image  pixel  with  the  highest  and  lowest  intensity  to  white  and 
black,  respectively.  The  camera  is  remotely  controlled  using  the  program  HyperTerminal 
from  Windows  and  its  focus  is  manually  adjusted  for  each  picture  until  the  best  setting  is 
achieved.  In  our  work,  all  images  were  digitized  using  eight  bits  per  pixel  and  contained 
one  face  only.  The  noise  equivalent  difference  temperature  (NEDT)  of  the  camera  is 
about  60  mK.  The  transmission  time  per  frame  using  the  RS-232  interface  is  about  6  sec¬ 
onds. 


2.  Image  Acquisition 

The  image  acquisition  scheme  is  illustrated  in  Figures  2  and  3.  There  were  nine 
numbered  points  marked  on  the  walls.  Numbers  vary  from  1  to  9,  as  depicted  in  Figure  3. 
Each  volunteer  was  asked  to  rotate  his  or  her  head  in  the  direction  of  the  numbers.  An 
additional  picture  was  taken  by  asking  volunteers  to  look  at  a  random  place  within  the 
square  formed  by  the  extreme  marked  points.  All  volunteers  were  adults  (13  males  and  1 
female),  all  males  were  clean-shaven  except  for  1  subject  with  a  mustache.  For  each  vol¬ 
unteer  three  sets  of  images  with  different  expressions  were  obtained.  Each  set  contained 
10  images,  resulting  in  420  images.  A  sample  of  the  images  obtained  is  shown  in  Figure 
4. 
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Figure  2.  Image  acquisition  (lateral  view). 


Figure  3.  Image  acquisition  (front  view). 
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Figure  4.  Sample  of  the  images  obtained. 


3.  Image  Nomenclature  Convention 

The  images  obtained  were  stored  in  .pgm  format  using  the  following  nomencla¬ 
ture:  xx-yy-zz.pgm,  where  xx  is  the  person  number,  yy  represents  the  visual  direction  that 
the  person  should  be  looking  at  and  zz  is  the  section  number. 

The  range  of  parameters  selected  for  this  study  is: 

•  [1-6,8,9,11-16]  for  xx, 

•  [1-10]  for  yy, 

•  [1,5,6]  for  zz. 

The  section  number  (zz)  was  defined  as  follows: 

1-  for  a  neutral  expression  and  no  glasses. 
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5-  for  a  “smiling”  expression  and  no  glasses. 

6-  for  a  vowel  “u”  pronouncing  expression  and  no  glasses. 

For  example,  2-5-5.pmg  refers  to  subject  number  2,  looking  at  visual  direction  5, 
with  a  “smiling”  expression  and  no  glasses. 

This  Chapter  presented  the  infrared  imaging  concepts,  the  setup  used  to  collect 
the  images  and  the  image  nomenclature  convention.  Next  chapter  discusses  the  classifica¬ 
tion  schemes  investigated. 
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IV.  FACE  RECOGNITION 


Automatic  face  recognition  is  a  procedure  for  associating  a  face  to  an  identity.  We 
may  say  that  a  face  is  recognized  when  its  image  matches  reasonably  well  the  face  image 
of  someone  whose  identity  we  know.  This  procedure  may  be  divided  in  three  sequential 
steps:  a)  Image  acquisition;  b)  Face  detection  and  segregation;  and  c)  Face  classification, 
as  shown  in  Figure  5. 


Figure  5.  Automatic  face  recognition  procedure  block  diagram. 


A.  FACE  DETECTION  AND  IMAGE  SEGMENTATION 

IR  images  were  collected  with  the  IR-160  camera,  as  described  in  Chapter  II. 
First,  the  images  were  cropped  to  isolate  face-only  portions  to  minimize  the  potential  im¬ 
pact  due  to  different  backgrounds.  Next,  the  resulting  cropped  images  were  used  in  the 
training  and  testing  phases  of  the  automatic  recognition  system.  Two  methods  for  crop¬ 
ping  were  used,  manual  cropping  and  automatic  cropping.  These  methods  are  described 
next. 


1.  Manual  Face  Detection  and  Segmentation 

In  this  work,  the  face  contained  in  each  picture  was  detected  by  visual  inspection 
and  manually  cropped  using  a  rectangle  of  60  vertical  by  45  horizontal  pixels.  The  size  of 
this  rectangle  was  obtained  by  trial  and  error  so  that  it  included  most  of  the  facial  features 
of  each  subject,  leaving  very  little  visible  background.  Faces  were  cropped  just  above  the 
chin  and  above  the  eyebrows.  Next,  each  cropped  image  was  saved  in  .bmp  format  using 
the  same  nomenclature  as  the  image  it  was  obtained  from,  except  that  an  “-a”  was  ap- 
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pended  to  the  file  name,  such  as  2-3-4-a.pgm  for  example.  A  sample  cropped  image  is 
shown  in  Figure  6  below. 


Figure  6.  Cropped  IR  image. 

2.  Automatic  Face  Detection  and  Segmentation 

An  automated  face  detection  and  cropping  routine  was  also  implemented  to  de¬ 
crease  the  manual  processing  time.  It  operated  in  two  steps: 

•  First,  images  are  vertically  cropped  using  a  basic  edge  detection  scheme. 

•  Second,  a  fine  search  on  the  resulting  image  is  conducted  to  extract  the 
face-only  portion  of  the  image. 

a.  Face  Edge  Location 

Face  edge  location  is  conducted  in  the  following  three  steps.  First,  a  “dif¬ 
ference”  image  is  generated  by  subtracting  a  shifted  copy  of  the  image  from  the  original 
image,  resulting  in  an  “edge  enhanced”  image,  as  shown  in  Figures  7  and  8  below.  Pixel 
values  obtained  for  row  60  of  the  difference  image  are  plotted  in  Figure  9,  which  shows 
that  vertical  left  and  right  edges  may  be  easily  located  by  extracting  the  indices  corre¬ 
sponding  to  the  maximum  (-40)  and  minimum  (-117)  values  respectively.  The  overall 
vertical  left  and  right  face  edges  are  established  by  selecting  the  smallest  left  edge  index 
and  the  largest  right  edge  index  out  of  three  rows  selected  in  the  image.  The  resulting 
cropped  image  is  the  original  image  cropped  on  the  right  and  left  edges  of  the  face,  as 
illustrated  in  Figure  10. 
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Figure  7.  Infrared  image. 


Figure  8.  Difference  image  obtained  from  Figure  7. 


Figure  9.  Row  number  60  pixel  values  obtained  from  Figure  8. 
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Figure  10.  Automatically  cropped  image. 

b.  Face  Locating  Procedure 

At  this  stage  we  still  need  to  further  crop  the  top  and  bottom  portions  of 
the  cropped  image  to  isolate  the  face-only  section  and  remove  any  left-over  background. 
Turk  and  Pentland  proposed  a  procedure  to  extract  exact  face  location  contained  in  an 
image  [30].  Their  algorithm  uses  the  Principal  Component  Analysis  (PC A)  approach, 
which  is  explained  in  further  details  in  the  next  section,  and  considers  the  Euclidian  dis¬ 
tance  as  a  measure  of  the  “closeness”  of  two  images.  Basically,  their  procedure  divides 
the  image  into  sub-images,  projects  them  using  PC  A  into  a  predefined  “face  space”  and 
uses  the  fact  that  “face-only”  sub-images  are  closer  to  the  face  space  than  non  “face- 
only”  subimages  are.  The  procedure  is  described  next. 

First,  we  extract  sub-images  r  from  an  image  /  and  subtract  the  mean 
image  of  the  whole  database  W,  as  the  mean  does  not  contain  any  useful  infonnation  for 
this  task.  The  resulting  mean  centered  sub-images  are  called  <p.  Next,  the  mean-centered 
sub-images  are  projected  onto  the  face  space,  which  leads  to  their  associated  sub-images 
<j)f.  At  that  point,  we  use  the  Euclidian  distance  to  evaluate  how  close  each  sub-image  is 
from  its  projection  onto  the  face  space  and  select  as  the  “face-only”  sub-image  that  which 
leads  to  the  minimum  distance  sr  out  of  all  sub-images  extracted  from  that  image.  Note 
that  this  method  is  computationally  expensive.  However,  we  can  obtain  the  face  map  di¬ 
rectly  from  the  original  image  without  explicitly  applying  the  projection  operation,  as 
shown  in  [30,  Eq.  15].  The  distance  expression  sr  can  be  shown  to  be  equal  to: 
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(x,  y)  =  r  (x,  y)T(x,  y)  -  2T(x,  y)  0 V  +  W  -  £  [r'  (x,  v)  0  m,  -  Y  0  u,  f,  (4.1) 

i= 1 

where  /"is  the  subimage,  0  is  the  correlation  operation,  W  is  the  mean  image  of  the 
whole  database  used  for  the  system  design  and  w,  are  the  eigenvectors  of  the  covariance 
matrix  used  to  create  the  face  space,  as  defined  in  the  PCA  section  later.  The  computation 
of  the  face  map  involves  only  the  computation  of  L+\  correlation  terms  over  the  input 
image,  and  the  computation  of  T‘(x,y)T(x,y) .  Computing  the  expression 
r'(x,y)r(x,y)is  obtained  by  squaring  the  input  image  I(x,y)  and  summing  the  squared 
values  contained  in  each  subimage.  This  procedure  was  implemented  in  MATLAB  and 
the  code  is  included  in  the  Appendix. 


B.  FACE  CLASSIFICATION 

Classification  aims  to  assign  an  identity  to  an  input  facial  image  by  comparing  it 
to  pre-stored  facial  images.  A  direct  image-by-image  comparison  would  be  extremely 
computationally  expensive.  Therefore,  a  dimension  reduction  process,  which  reduces  the 
search  time  and  memory  required  to  store  all  subjects,  is  very  desirable.  Note  that  an  im¬ 
age  may  be  considered  as  a  vector  of  dimension  equal  to  the  number  of  pixels  in  the  im¬ 
age  by  reshaping  the  image  as  a  vector.  For  example,  consider  a  60  by  45  pixels  image, 
which  contains  a  total  of  2,700  pixels.  As  a  result,  the  image  vector  produced  by  this  im¬ 
age  is  of  length  2,700.  In  addition,  pixels  are  likely  to  be  correlated.  As  a  result,  a  dimen¬ 
sion  reduction  scheme  is  likely  to  reduce  computational  load  and  storage  requirements 
without  degrading  discrimination  performances.  In  this  work,  we  consider  two  linear  pro¬ 
jection  approaches  which  allow  us  to  significantly  reduce  the  problem  dimension:  the  Ei- 
genface  approach,  based  on  the  Principal  Component  Analysis  (PCA)  concept,  and  the 
Fisherface  approach,  based  on  the  PCA  and  Fisher  Linear  Discriminant  (FLD)  concepts. 
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1.  Principal  Component  Analysis  (PCA) 
a.  Definition 

PCA  is  a  linear  projection  scheme  that  reduces  the  dimensionality  of  a 
data  set  while  retaining  as  much  variance  present  in  the  data  as  possible.  As  a  result,  the 
PCA  scheme  is  best  matched  for  compression  applications,  but  has  also  been  applied  with 
success  to  classification  problems,  although  one  cannot  insure  the  components  kept  are 
those  best  suited  for  discrimination  purposes  [31]. 

Initially,  let  us  consider  the  problem  of  representing  N  ^-dimensional 
samples,  image  vectors  {xj,  xj,  ...,  x a?}  by  a  single  vector  xq.  To  be  more  specific,  sup¬ 
pose  that  we  want  to  find  a  vector  xo  such  that  the  sum  of  the  squared  distances  between 
x_o  and  the  various  xj;,  k  =  is  as  small  as  possible.  Let  us  define  a  cost  function 

Jo(xo)  that  must  be  minimized  as: 

n 

*4(Lo)  =  Zfe>-&|f  •  (4-2) 

k= 1 


Next,  we  will  show  that  the  solution  to  this  problem  is  given  by  x_o=m, 
where  m  is  the  dataset  sample  mean, 


1  n 

Ul  =  ~Y*k  • 
n  /,=, 


(4.3) 


First,  let  us  rewrite  Jo(xo)  as  follows, 

J0(x o)  =  X|(Lo  - m)-(xk  -m)\\ 

k=l 

n  ^  n  n 

=  X  It  -  -  2X  (To  -  m)'  (**  -m)  +  Yj  It  _  1 1 

k= 1  k= 1  k= 1 

n  n  n 

=  Yj  It  -  ™(  -  2(l0  ~m)'Y  (**  ~m)+Y  It _  1 1 

k= 1  k= 1 

n  ,  ^  w  ^  2 

=it-tr+iit-tr- 


k= 1 


(4.4) 


k= 1 


k= 1 
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Note  that  the  above  expression  is  minimized  when  xo=m,  since  the  second 
tenn  in  the  summation  is  independent  of  x_o. 

Next,  let  us  consider  e  to  be  a  unit  vector  in  the  direction  of  the  line  run¬ 
ning  through  the  sample  mean.  The  equation  of  the  line  can  be  written  as: 

x  =  m  +  ae,  (4.5) 


where  the  scalar  a  (which  can  take  any  real  value)  corresponds  to  the  distance  of  any 
point  x  from  the  mean  m.  At  this  stage,  we  can  find  an  “optimal”  set  of  coefficients  a*’ s, 
if  we  represent  xy  by  m+ai<e,  by  minimizing  the  following  cost  function: 


Jx  =  (a1,a2,...,a„,e)  =  ^|(m  +  aite)-xt 

k= 1 


Xlh :e-(Xk~m)\ 

k= 1 


=  Z 1 ak2  IMf  -  2Z  ak  l  (**  -m)  +  Y\xk-m\ 

k= 1  k= 1  k= 1 


(4.6) 


Recognizing  that  ||e||  is  equal  to  1,  partially  differentiating  with  respect  to  a*  and  setting 
the  derivative  to  zero,  leads  to: 

ak  =e{xk-m).  (4.7) 

Geometrically,  this  result  means  that  least-square  solution  may  be  obtained  by  projecting 
the  vector  xy  on  the  line  in  the  direction  of  e  that  passes  through  the  sample  mean  m. 

Next  let  us  consider  the  problem  of  finding  the  best  direction  e  for  the  line 
of  interest.  Let  us  define  the  total  scatter  matrix  S  as: 


S  =  ^d(xk-m)(xk-!n)t  ■ 

k= 1 

Replacing  the  Equation  (4.7)  into  Equation  (4.6)  leads  to: 

n  n  n 

=  Yjak2  ~2Yjak2  +  Zb- 

k= 1  k= 1  k= 1 

n  r  i2  11 


(4.8) 


k= 1 


k= 1 
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=  -JV  (**  “  *  -  «)'  e  +  X  fc*  -  5»|| 

Ar=l  £=1 

=  -g^g  +  yjb  -m||2.  (4.9) 

t=i 

Naturally,  the  vector  e  that  minimizes  Ji  also  maximizes  eSe.  The  solution 
to  this  problem  must  satisfy  the  equation: 

Se  =  Ae.  (4.10) 

This  means  that  the  vector  e  that  minimizes  Ji(e)  is  an  eigenvector  of  the  scatter  matrix. 
As  a  result,  maximizing  the  expression  eSe  is  obtained  by  selecting  for  e  the  eigenvector 
of  S  associated  with  the  largest  eigenvalue.  In  other  words,  we  project  the  data  onto  a  line 
running  through  the  sample  mean  in  the  direction  of  the  eigenvector  of  the  scatter  matrix 
having  the  largest  eigenvalue  to  find  the  best  one-dimensional  projection  of  the  data. 

Figure  1 1  presents  a  two-class,  two-dimensional  problem  where  we  wish 
to  reduce  the  data  dimension  from  two  to  one.  Class  1  and  2  are  represented  by  stars  and 
circles  respectively.  The  direction  that  provides  the  best  one-dimensional  PCA  represen¬ 
tation  of  the  data  is  the  line  included  in  Figurell.  Note  that  class  discrimination  is  pre¬ 
served  as  the  projected  stars  and  circles  onto  the  line  do  not  overlap. 
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Figure  1 1 .  PCA  Example.  The  line  represents  the  projection  direction  in  which  the  data 
should  be  projected  to  obtain  the  best  representation  in  one  dimension. 

Recall  that  PCA  was  initially  designed  for  compression  applications,  and 
that  the  projection  directions  are  not  optimized  to  preserve  class  separability.  Figure  12 
shows  an  example  in  which  the  projection  direction  that  best  represents  the  data  does  not 
preserve  class  discrimination  when  projecting  a  two-dimensional  dataset  with  two  classes 
onto  a  one-dimensional  space  (a  line).  In  this  example,  the  projected  stars  and  circles 
overlap,  and  it  is  no  longer  possible  to  separate  the  two  classes. 
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Figure  12.  Data  Points  and  PC  A  projection.  The  line  provides  the  best  projection  direction  to 

represent  the  data  set. 

The  above  result  can  be  extended  from  a  one-dimensional  projection  to  a 
/^-dimensional  projection. 

p 

x  =  tn  +  '5',ai g,:  •  (4-11) 

1=1 


The  cost  function  becomes 


n 

I 

k= 1 


v 


(4.12) 


This  function  is  minimized  when  the  vectors  e_i,  e2,...ep  are  the  p  scatter 
matrix  eigenvectors  associated  with  the  p  largest  eigenvalues.  Note  that  these  eigenvec¬ 
tors  are  orthonormal  as  the  scatter  matrix  is  real  and  symmetric,  and  they  form  a  basis  set 
used  to  represent  any  feature  vector  x.  The  coefficients  cik  are  the  projections  of  x  in  that 
basis  and  are  called  principal  components. 
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The  PC  A  scheme  may  be  implemented  in  a  matrix  fonn,  as  follows: 

Each  image  contained  in  the  training  set  is  reshaped  as  a  column  vector 
with  dimension  q,  and  the  data  matrix  A  is  defined  as  the  concatenation  of  all  image  vec¬ 
tors  columnwise,  resulting  in  a  matrix  of  dimension  q  (in  our  case  2700).  Next,  the  mean 
image  of  the  training  set  is  subtracted  from  each  image  resulting  in  a  matrix  X.  The  data 
covariance  matrix  S  is  defined  as: 

S=XXt.  (4.13) 

Recall  that  the  d  eigenvectors  of  S  corresponding  to  the  largest  eigenval¬ 
ues  minimize  the  cost  function  ./(■)  previously  defined.  Further,  note  that  the  matrix  S 
has  dimension  2700  in  our  case,  which  is  very  large.  Computing  the  eigenvectors  of  S  = 
XX 4  may  be  significantly  decreased  by  computing  them  in  terms  of  the  eigenvectors  of 
XtX with  the  snapshot  method  which  is  described  next. 

b.  Snapshot  Method 

The  singular  value  decomposition  theorem  states  that  any  K  x  N  matrix 
can  be  decomposed  and  written  as  the  product  of  three  matrices  [32] 

X=  UZV*T ,  (4.14) 

where  U  is  the  K  x  K  unitary  left  singular  vectors  matrix,  containing  the  singular  vectors 
columnwise,  V  is  the  N  x  N  unitary  right  singular  vectors  matrix,  and  E  is  the  K  x  /V  di¬ 
agonal  matrix  of  nonnegative  real  singular  values.  The  correlation  matrix  XX*  can  be 
expressed  as: 

XX1  =  (UZV')(VZTU')=  UL{ytV)'L,Ut  =  ITLHU* .  (4.15) 

By  analogy  the  product  XtX can  be  expressed  as: 

X‘X=  VYlYV‘ .  (4.16) 

Equation  (4.15)  shows  that  the  matrix  U  contains  the  eigenvectors  of  XX1 

and  the  eigenvalues  of  XXt  are  diagonal  elements  of  the  matrix  EE' ,  by  analogy  Equation 
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(4.16)  shows  that  the  matrix  V  contains  the  eigenvectors  of  XX  and  that  the  eigenvalues 
of  XX  are  diagonal  elements  of  the  matrix  Z'Z . 

Next,  let  us  consider  the  product  XV, 

XV  =  UZV'V  =  UZ,  (4.17) 

where  the  last  step  follows  because  V  is  unitary.  Note  that  the  only  useful  eigenvectors 
for  our  application  are  those  associated  with  non-zero  eigenvalues.  Therefore  we  define 
I  1  as  the  inverse  of  the  portion  of  I  which  contains  the  non-zero  diagonal  elements  of 
Z  only.  Therefore,  the  eigenvectors  associated  with  non-zero  eigenvalues  are  defined: 

U=  XVZ'1.  (4.18) 


2.  Linear  Discriminant  Analysis  (LDA) 

LDA  is  an  approach  designed  to  reduce  a  dataset  dimensionality  while  maximiz¬ 
ing  class  separation.  As  a  result,  the  LDA  scheme  is  optimized  for  classification  prob¬ 
lems,  where  the  goal  is  to  preserve  class  discriminative  information  [31]. 

Let  us  consider  initially  the  case  of  finding  one  line  that  best  discriminates  the 
projections  of  a  two-class  J-dimensional  dataset  and  later  we  will  generalize  the  equation 
to  a  c  classes  problem. 

Assume  that  we  have  a  set  of  n  ^-dimensional  samples  xj,  x2,  ...xn,  where  n /  sam¬ 
ples  are  in  subset  (i.e.,  class)  C/,  and  n2  samples  are  in  subset  (i.e.,  class)  C2.  If  we  form  a 
linear  combination  of  the  components  of  x  we  obtain: 

y  =  wx,  (4.19) 

and  a  corresponding  setyy,  y2,...,y„  divided  in  two  subsets  Y/  and  Y2. 

Geometrically,  if  ||w||  is  equal  to  one,  each  v,  is  the  projection  of  the  correspond¬ 
ing  Xj  onto  a  line  in  the  direction  of  w.  Note  that  the  magnitude  of  w  is  of  no  real  signifi¬ 
cance  because  it  merely  scales  y.  A  measure  of  the  separation  between  the  projected 
points  is  the  difference  of  the  sample  means. 
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Let  us  define  m,  the  ^/-dimensional  sample  mean  as  follows: 


"L  =  — Z*’ 


(4.20) 


XGX; 


where  n,  is  the  number  of  elements  on  the  set  X,. 

The  sample  mean  for  the  projected  points  is  given  by 


q,  = 


y*=Y, 


—  7  w  x  =  w  mt. 

ni  ieD. 


(4.21) 


Note  that  q  (above)  is  simply  the  projection  of  m,on  w. 


The  distance  between  the  projected  means  is  given  by: 


-<h\ 


w(mx  -m2) 


(4.22) 


The  above  distance  may  be  varied  by  scaling  the  vector  w.  Thus,  a  good  measure 
of  separation  for  the  projected  data  requires  the  definition  of  this  distance  with  respect  to 
some  measure  of  the  standard  deviations  for  each  class. 


We  can  define  the  scatter  for  projected  samples  by: 

(A)2  =  2>-?,)2  ■  (4.23) 

y&Yi 


The  Fisher  Linear  discriminant  employs  the  linear  transformation  w  such  that  the  cost 
function: 


J(w)  = 


m  -<h\ 


(a)2 -(4  r 


(4.24) 


is  maximized.  The  direction  w  maximizing  ./(•)  leads  to  the  best  separation  between  the 
two  projected  sets.  To  obtain  ./(■)  as  a  function  of  w,  we  define  the  scatter  matrices  S,- 
and  Sw  by: 

S;  =  Z  y  ,  (4.25) 


and 
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(4.26) 


SW=S1+S2. 


Rewriting  Lt  from  Equation  (4.23),  the  scatter  from  the  projected  data,  we  have: 

(A  )  =  X!  Vm,)2  =  y'i  w  (x - mt )(x -m,yw 

x&Xj  xeXj 

=  wS,w.  (4.27) 

Therefore  the  sum  of  these  scatter  matrices  can  be  written 

(f1)2+(i2)2=4w.  (4.28) 

Similarly,  the  separations  of  the  projected  means  is  given  by: 

(g,  ~q2)  =  (w*Z?h -wtJh) 

=  w‘SBw,  (4.29) 

where  SB  =(ml  - m2)(m ,  -m2)‘ .  (4.30) 

Sw  is  the  within-class  scatter  matrix.  It  is  symmetric  and  positive  definite,  and  it  is 
usually  nonsingular  if  n>d.  Likewise,  S/S  is  called  the  between-class  scatter  matrix.  It  is 
also  symmetric  and  positive  definite,  and  has  rank  one  because  it  is  defined  as  the  outer 
product  of  two  vectors.  In  particular,  for  any  w,  Si,w  lies  in  the  direction  of  m_i  -mj. 

The  cost  function  J(-)  can  be  expressed  in  tenns  of  Sb  and  Sw  as: 


J(w)  = 


wSBw 

wSww 


(4.31) 


Therefore,  the  vector  w  that  maximizes  ./(■)  must  satisfy  SBw  =  ASww,  which  is  a  gen¬ 
eralized  eigenvalues  problem.  This  generalized  eigenvalue  problem  may  be  reformulated 
as  conventional  eigenvalue  problem  when  Sw  is  nonsingular,  writing 

S;;SHw  =  Xw.  (4.32) 
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In  our  particular  case,  it  is  unnecessary  to  solve  for  the  eigenvalues  and  eigenvec¬ 
tors  of  S~1SB  due  to  the  fact  that  SBw  is  always  in  the  direction  of  mi-mj.  As  a  result, 


vector  w  that  optimizes  /(•)  is  given  by: 

w  =  Sn'(ml  -m2). 


(4.33) 


As  an  example  consider  the  two-dimensional  two-class  data  points  shown  in  Fig¬ 
ure  13,  where  we  are  interested  in  the  one-dimensional  projection  which  preserves  the 
class  separation.  The  vector  w  that  maximizes  the  discriminant  between  the  two  classes 
(stars  and  circles)  is  the  vertical  line  which  shows  the  projected  data  can  still  be  discrimi¬ 
nated  on  the  line. 
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Figure  13.  Data  points  and  LDA.  The  line  provides  the  best  projection  direction  to  discrimi¬ 
nate  the  data  sets. 
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The  linear  discriminant  can  be  extended  to  a  c-class  problem.  The  natural  gener¬ 
alization  involves  c-1  discriminant  functions.  Thus,  the  projection  is  from  a  d- 
dimensional  space  to  a  (c-  1)  dimensional  space,  where  it  is  assumed  that  d  >  c  . 

The  within-class  scatter  matrix  for  the  c-class  problem  is  defined  as: 

Sw=fJSi,  (4.34) 

i= 1 

where, 

S,  =  X  (*“ (4.35) 

xeXi 

and 

mj  =  y'y  x  .  (4.36) 

xeXj 

The  proper  generalization  for  SB  is  not  quite  so  obvious.  Suppose  that  we  define  a 
total  mean  vector  m  and  a  total  scatter  matrix  St  by: 

m  =  —  ^  x,  (4.37) 

n  x  g  all  xt 

and 

ST  =y'Xx~ m)(x  - m)'  .  (4.38) 
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The  projection  from  a  t/-dimensional  space  to  a  (c-1)  dimensional  space  is  ac¬ 
complished  by  c - 1  discriminant  functions  as:  yi=wjx  ,  i  =  1,...,  c-1.  Note  that  the 

parameters  y,  may  be  viewed  as  the  components  of  a  vector  y  and  the  weight  vectors  w, 
may  be  viewed  as  the  columns  of  a  d-by-c  matrix  W.  In  such  a  case  the  projection  opera¬ 
tion  can  be  rewritten  as  a  single  matrix  equation: 

y  =  W‘x.  (4.40) 

Thus,  the  samples  xi,  ...  ,  xn  project  to  a  corresponding  set  of  samples  y/,  ...,y, 
which  can  be  described  by  their  own  mean  vectors  and  scatter  matrices. 

Let  us  define. 


and 


It  can  be  proved  [31]  that 


and 


q,  =-I >, 

>h  y&t  ~ 

(4.41) 

q  =  -  Z  y  ’ 

n  ysallYi 

(4.42) 

:SXcv-9,.)(v-?.y, 

i=l  ytY; 

c 

(4.43) 

YJnMi  -q  )(q,  -q  )'  • 

i= 1 

(4.44) 

Lw  =  W‘SwW , 

(4.45) 

lb  =  w'sBw  . 

(4.46) 

At  this  point,  we  seek  is  a  transformation  matrix  W  that  in  some  sense  maximizes 
the  ratio  of  the  between-class  scatter  to  the  within-class  scatter.  A  simple  scalar  measure 
of  the  scatter  is  the  detenninant  of  the  scatter  matrix.  Recall  that  the  determinant  is  the 
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product  of  the  eigenvalues  and,  hence,  is  the  product  of  the  “variances”  in  the  principal 
direction. 


Using  this  measure,  we  obtain  the  cost  function: 


J(W)  = 


Lb 

W‘SBW 

Lw 

W'SwW 

(4.47) 


The  columns  of  an  optimal  W,  i.e.,  defined  as  maximizing  the  cost  function  ,/(■) , 
corresponds  to  the  largest  eigenvalues  of  the  generalized  eigenproblem: 


SBwt  =  A,Sww,. . 


(4.48) 


Note  that  the  within-class  scatter  matrix  SV  is  always  singular  in  the  face  recogni¬ 
tion  problem.  This  singularity  constraint  is  due  to  the  fact  that  the  rank  of  Sw  is  at 
most  N  -c  ,  where  N  is  the  number  of  images,  c  is  the  number  of  classes  and  the  number 
of  images  is  smaller  than  the  number  of  pixel  in  each  image  n. 

In  order  to  overcome  this  complication  of  a  singular  Sw,  an  alternate  cost  function 
proposed  by  Belhumeur,  Hespanha  and  Kriegman  [33]  was  used  in  this  work.  This 
method,  called  Fisherface,  avoids  the  singularity  problem  by  projecting  the  image  set  to  a 
lower  dimensional  space  so  that  the  resulting  within-class  scatter  matrix  Sw  is  nonsingu¬ 
lar.  The  dimension  reduction  is  achieved  by  using  PCA  to  reduce  the  feature  space  di¬ 
mension  to  N  -c ,  where  N  it  the  total  number  of  images  used  on  the  training  set,  and 
then  applying  the  standard  LDA  to  further  reduce  the  dimension  to  c  - 1 . 


3.  Minimum  Distance  Classifier 

Once  projection  directions  are  computed  from  the  training  set,  class-specific  cen¬ 
troids  are  obtained  by  computing  the  average  values  of  the  projected  data  for  each  class. 
Next  the  testing  images  are  projected  using  the  projection  directions  defined  for  the  train¬ 
ing  set,  and  the  distance  between  this  projected  vector  and  each  class  centroids  is  com¬ 
puted.  The  classification  decision  is  obtained  by  selecting  the  class  that  is  the  closest  in 
nonn  to  the  projected  testing  image. 
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This  Chapter  presented  the  two  linear  dimension  reduction  schemes  implemented, 
the  snapshot  method  implemented  to  reduce  the  computational  cost  associated  with  the 
PCA-based  algorithms,  and  described  the  minimum  distance  classifier  implemented  in 
this  study.  Next  Chapter  presents  the  classification  results  obtained  with  our  IR  database. 
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V.  RESULTS 


PCA  and  LDA-based  face  recognition  system  classification  performances  are  es¬ 
timated  using  a  cross-validation  variant  to  take  into  account  the  relatively  small  size  of 
the  database  (420  images  corresponding  to  14  individuals  with  30  images  each  were  gen¬ 
erated  and  used  in  the  results). 


A.  CROSS  VALIDATION 

Cross-validation  is  a  statistical  scheme  designed  to  estimate  generalization  errors 
based  on  "resampling".  It  is  quite  useful  when  dealing  with  small  databases  as  it  aver¬ 
ages  classification  performances  over  multiple  testing  subsets  [34].  The  basic  k-fold 
cross-validation  implementation  divides  the  dataset  into  k  disjoint  sets,  where  training  is 
conducted  with  k  - 1  datasets  and  tested  with  the  remaining  set.  The  process  is  repeated  k 
times,  where  each  time  the  selected  testing  set  is  different,  and  the  overall  classification 
performance  is  obtained  as  the  mean  of  all  perfonnances  found. 

In  this  work,  we  consider  a  cross-validation  variant,  where  the  data  available  is 
separated  in  two  disjoint  sets,  the  training  and  the  validation  set,  as  shown  in  Figure  14. 
The  training  set  containing  60%  of  the  pictures  available  is  used  to  obtain  the  projection 
directions  and  class  centroids.  The  validation  set  contains  the  remaining  40%  of  the  pic¬ 
tures  available  and  is  used  to  test  the  system  performance.  This  data  split  and  classifica¬ 
tion  estimation  process  is  repeated  1000  times,  while  the  number  of  images  per  subject  is 
kept  constant  for  all  subjects,  and  the  overall  classification  performance  is  computed  as 
the  mean  of  all  repetitions. 
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Figure  14.  Block  diagrams  illustrating  the  cross  validation  procedure. 

B.  PCA  RESULTS 

In  addition  to  the  basic  PCA-based  implementations,  we  also  considered  the  im¬ 
pact  of  removing  the  first  few  eigenvector  contributions,  as  is  done  successfully  in  visible 
imagery  classification  to  make  the  recognition  algorithms  more  robust  to  illumination 
changes  [33], 

Several  PCA-based  variants  were  investigated: 

•  PCA  100:  PCA  scheme  using  the  first  top  100  eigenvectors  for  projection 
directions  (i.e.,  the  eigenvectors  associated  with  the  100  largest  eigenval¬ 
ues); 

•  PCA40:  PCA  scheme  using  the  first  top  40  eigenvectors  for  projection  di¬ 
rections; 

•  PCA40W1:  PCA  scheme  using  the  first  40  eigenvectors,  after  removing 
the  top  first  eigenvector; 

•  PCA40W2:  PCA  scheme  using  the  first  40  eigenvectors,  after  removing 
the  top  two  eigenvectors; 

•  PCA40W3:  PCA  scheme  using  the  first  40  eigenvectors,  after  removing 
the  top  three  eigenvectors; 

•  PCA40W4:  PCA  scheme  using  the  first  40  eigenvectors,  after  removing 
the  top  four  eigenvectors; 
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•  PCA40W5:  PCA  scheme  using  the  first  40  eigenvectors,  after  removing 
the  top  five  eigenvectors. 

The  maximum  number  of  PCA-based  projection  directions  is  related  to  the  di¬ 
mension  of  the  data  used,  but  some  projection  directions  may  be  disregarded  because  the 
corresponding  eigenvalues  associated  are  close  to  zero.  It  can  be  shown  that  the  maxi¬ 
mum  number  of  PCA-based  projection  directions  associated  with  non-zero  eigenvalues  is 
equal  to  the  size  of  the  dataset  when  the  database  size  is  smaller  than  the  data  dimension. 
In  our  work  the  maximum  number  of  PCA-based  projection  directions  is  2700,  and  the 
number  of  PCA-based  projection  directions  associated  with  non-zero  eigenvalues  is  256 
(that  is,  equal  to  60%  of  420).  The  number  of  “useful”  PCA-based  projection  direction 
was  found  experimentally  to  be  equal  100.  Note  that  we  defined  a  projection  direction 
(i.e.,  eigenvector)  as  “useful”  when  the  associated  eigenvalue  is  at  least  equal  to  0.1%  of 
the  largest  value.  Figure  15  plots  the  overall  classification  performance  when  a  combina¬ 
tion  of  datasets  #1  and  #6  were  used  as  the  training  set,  and  dataset  #5  was  selected  for 
testing.  Figure  15  shows  that  the  best  classification  performance  that  can  be  achieved 
with  the  data  used  in  this  study  is  obtained  with  the  top  thirty-one  eigenvectors  (leading 
to  an  error  rate  of  26.5%).  Note  however,  that  this  result  is  for  one  training/testing  split 
combination.  Similar  behavior  was  observed  for  other  combinations.  Therefore,  we  se¬ 
lected  the  top  40  eigenvectors  to  insure  the  best  obtainable  classification  performances 
for  all  training/testing  sets  combinations  used  for  the  cross-validation  implementation. 

The  cross-validation  implementation  uses  1000  repetitions,  where  the  datasets  are 
kept  identical  for  all  methods  used  in  each  case.  Classification  error  histograms  for  all 
schemes  are  shown  in  Figures  16  to  24.  Figure  25  illustrates  the  LDA  classification 
scheme  robustness,  as  explained  further  below.  Table  2  lists  the  mean  classification  error 
rate  in  percent  for  each  subject  considering  the  1000  repetitions. 
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Figure  15.  Error  rate  (%)  versus  the  number  of  eigenvectors 


Figure  16.  LDA  error  rate  histogram  in  %;  cross-validation  implemented 
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%  of  Simulations  erg’  %  of  Simulations 


are  17.  Direct  classification  error  rate  histogram  in  %;  cross-validation  implemented. 
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%  of  Simulations 


Figure  19.  PCA40  error  rate  histogram  in  %;  cross  validation  implemented. 


tion  implemented. 
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Figure  2 1 .  PCA40W2  error  rate  histogram  in  %  with  first  two  eigenvectors  removed;  cross 

validation  implemented. 


Figure  22.  PCA40W3  error  rate  histogram  in  %  with  first  three  eigenvector  removed;  cross 

validation  implemented. 
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%  of  Simulations 


Figure  23.  PCA40W4  error  rate  histogram  in  %  with  first  four  eigenvector  removed;  cross 

validation  implemented. 
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Figure  24.  PCA40W5  error  rate  histogram  in  %  with  the  first  five  eigenvectors  removed; 

cross  validation  implemented. 


Person 

Figure  25.  LDA  scheme;  Minimum/Mean/Maximum  class  distances  from  testing  set  subject 
#1  to  all  training  set  class-specific  centroids;  training  set:  datasets  #1  and  #5;  test¬ 
ing  set:  dataset  #6. 
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Table  2.  Error  rate  (%)  per  class. 


Class 

1 

2 

3 

4 

5 

6 

8 

9 

11 

12 

13 

14 

15 

16 

Mean 

Error 

Rate  (%) 

Method 

Direct 

42 . 62 

36.08 

11.08 

14 . 07 

0.01 

16.01 

11.93 

0.01 

10.05 

3.34 

75 . 72 

3.29 

0.19 

2.76 

16.22 

PCA100 

42 . 65 

36.09 

11.18 

14  .  17 

0.01 

16.13 

11.97 

0.01 

10.15 

3.34 

75.79 

3.29 

0.20 

2 .80 

16.27 

PCA4  0 

42 . 97 

36.25 

11.93 

14 . 75 

0.01 

16.63 

12.10 

0.01 

10.99 

3.38 

75.88 

3.33 

0.30 

3.58 

16.58 

PCA40W1 

19.87 

7 . 93 

26.23 

0.00 

30.37 

11.39 

22 .45 

0.00 

8 . 68 

3.93 

30.49 

0 . 08 

9 . 65 

5.02 

12.58 

PCA40W2 

15 . 12 

11.24 

13.03 

0.00 

30.03 

12 . 85 

23.06 

0.04 

10.53 

3.32 

33.84 

0.76 

6.39 

3.69 

11 . 71 

PCA40W3 

15.11 

9.06 

7 . 51 

0.00 

17 .40 

3.21 

15 . 98 

0 . 04 

5.11 

3 . 93 

24.31 

0.26 

1 . 99 

5 . 08 

7 . 78 

PCA40W4 

22 . 18 

10 . 00 

22 . 52 

0.00 

11.40 

5 . 08 

14.93 

3.11 

3 . 92 

4.78 

32.56 

0.54 

0 . 05 

11.13 

10 .16 

PCA40W5 

23.52 

12 . 72 

33.63 

0.07 

22.39 

8 . 60 

20.59 

21.48 

6.30 

10.57 

21.21 

3 . 08 

0 .48 

13 .16 

14 . 13 

LDA 

2 . 50 

0 . 00 

1.91 

0.00 

0.03 

1 .75 

0.30 

0 . 02 

0.50 

0 . 00 

0 . 83 

0 . 00 

0 . 12 

0.50 

0 . 60 

The  following  comments  can  be  made: 

1.  Figure  16  shows  the  results  for  the  LDA  scheme,  it  provides  an  overall  average 
classification  error  rate  equal  to  0.62%.  On  36%  of  the  repetitions  there  were 
no  classifications  errors  and  on  35%  of  the  simulations  there  was  1  classifica¬ 
tion  error. 

2.  Figure  17  shows  the  results  for  the  direct  classification  without  the  use  of  di¬ 
mension  reduction.  On  Figures  18  and  19  are  depicted  the  results  for  the 
PCA100  and  PCA40  schemes.  From  these  figures  we  verify  that  there  were  no 
degradation  on  the  system  performance  as  we  reduce  the  dimension  of  the  data 
before  applying  the  minimal  distance  classifier. 

3.  From  Figures  20,  21  and  22  we  verify  that  the  mean  error  rate  decreases  as  we 
remove  the  first  top  three  eigenvectors.  From  Figures  23  and  24  we  verify  that 
the  mean  error  rate  increases  if  we  eliminate  further  top  eigenvectors.  The  best 
PCA-based  scheme  was  the  PCAW3  that  has  an  overall  classification  error  rate 
equal  to  7.78%. 

4.  Figure  25  illustrates  the  LDA-based  classification  perfonnance  robustness  by 
plotting  the  distances  between  one  image  of  test  subject  #1  to  all  16  training  set 
centroids,  when  datasets  #1  and  #5  are  selected  for  training  and  dataset  #6  is 
used  for  testing.  The  bar  plot  shows  minimum,  maximum  and  average  dis¬ 
tances  from  test  subject  #1  to  any  class-specific  centroid.  Note  that  the  dis¬ 
tance  from  test  subject  #1  to  the  class  #1  centroid  is  significantly  smaller  than 
it  is  to  any  other  class  centroid,  indicating  the  class-specific  clusters  are  well 
separated  as  far  as  the  test  subject  is  considered,  which  makes  the  classification 
decision  easy  to  make.  Similar  behavior  can  be  observed  for  all  other  test  sub¬ 
jects  in  this  testing  set. 

5.  The  LDA-based  classification  performance  was  evaluated  by  removing  the 
largest  eigenvector  but  this  step  did  not  result  in  any  improvement. 

6.  Table  2  shows  the  performance  of  each  classification  method  for  each  set  of 
images  (i.e.,  each  “class”  of  images).  The  results  show  the  sensitivity  of  each 
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class  of  images  to  the  classification  method.  For  example,  some  classes  showed 
improvement  with  the  removal  of  the  top  eigenvector;  others  did  not.  From  the 
right  column  of  the  table,  we  see  that  the  LDA  method  had  the  lowest  mean  er¬ 
ror  rate.  The  PCA40W3  had  the  lowest  mean  error  rate  of  all  of  the  PCA-based 
methods. 


This  Chapter  presented  the  face  recognition  classification  performances  obtained 
with  the  two  dimension  reduction  approaches  considered  in  this  study.  Results  show  that 
best  performances  are  obtained  with  the  Linear  Discriminant  implementation,  as  expected 
from  its  definition  that  is  best  matched  for  classification.  Results  also  show  that  the 
PCA-based  implementation  benefits  from  removing  the  top  three  eigenvectors,  pointing 
out  the  fact  that  there  are  still  brightness  issues  which  need  to  be  addressed  when  dealing 
with  IR  imagery. 
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VI.  CONCLUSIONS 


This  study  investigated  a  (IR)  face  recognition  system  using  an  uncooled  IR  cam¬ 
era.  A  computer-based  image  collection  set-up  was  designed  and  used  to  collect  a  small 
database.  The  small  database  has  420  facial  images,  from  14  volunteers  (13  adult  males 
and  1  adult  female)  with  30  pictures  per  person.  All  subjects  are  without  glasses.  Three 
different  facial  expression  sets  were  collected;  first,  the  individuals  wore  a  neutral  ex¬ 
pression.  The  second  set  was  constituted  of  the  same  subjects  with  a  “smiling”  expres¬ 
sion.  The  third  set  was  collected  with  individuals  pronouncing  the  vowel  “u”.  Study  vari¬ 
ables  were  reduced  by  using  a  controlled  environment  with  the  same  background  for  each 
picture,  with  the  person-camera  distance  fixed,  and  by  restricting  the  pictures  to  frontal 
facial  images  while  allowing  a  vertical  and  horizontal  angle  freedom  of  10°. 

Manual  and  automated  facial  image  cropping  routines  were  implemented,  as  de¬ 
scribed  earlier  in  Chapter  IV.  The  cropping  format  was  a  simple  rectangle  of  fixed  size 
equal  to  60  by  45  pixels.  Two  linear  approaches  for  the  dataset  dimension  reduction  and 
classification  were  implemented  and  their  resulting  classification  perfonnances  com¬ 
pared,  PCA-based  and  LDA  approaches.  A  minimum  distance  classifier  was  selected  to 
evaluate  the  classification  performances.  The  overall  system  performance  was  evaluated 
with  a  cross-validation  scheme  using  60%  of  the  pictures  for  training  and  40%  of  the  pic¬ 
tures  for  testing,  with  1,000  repetitions. 

Results  show  that  the  best  PCA-based  overall  classification  perfonnance 
(92.22%)  is  obtained  when  selecting  the  top  40  eigenvectors,  while  excluding  the  first 
three  top  eigenvectors.  The  LDA-based  approach  perfonned  better,  with  an  overall  classi¬ 
fication  perfonnance  equal  to  99.40%,  as  expected  from  the  scheme  definition. 

Results  obtained  in  this  study  successfully  demonstrate  that  an  uncooled  IR  cam¬ 
era  may  discriminate  between  individual  subjects  obtained  from  a  small  database  col¬ 
lected  under  a  very  controlled  environment.  The  database  needs  to  be  extended  and  more 
sophisticated  image  normalization  investigated  to  take  into  account  facial  images  taken  at 
different  distances  and/or  resolutions.  In  addition,  this  study  considered  two  linear  pro- 
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jections  only.  More  sophisticated  nonlinear  implementations  are  expected  to  improve  the 
classification  performance  when  dealing  with  a  larger  database. 
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APPENDIX 


This  appendix  contains  all  MATLAB-based  software  implemented  for  this  study. 

Images  may  be  loaded  into  MATLAB  individually  or  in  batch  mode: 

•  Individual  loading  is  done  with  the  routine:  readpgm8.m. 

•  Batch  loading  is  done  with  the  routine:  Loadlmage.m,  where  the  operator 
selects  the  range  of  images  to  be  loaded  in  batch  mode. 

Images  may  be  loaded,  cropped  and  saved  back  into  .bmp  format  with  user- 
specified  names  in  batch  mode  using  the  function: 

•  CropFaceAuto.m  (calls  the  automated  face  location  code  described  be¬ 
low). 

Automated  face  location  extraction  code  is  implemented  with  the  functions: 

•  locatehead.m  which  vertically  crops  the  image  around  the  head, 

•  facemapini.m  which  computes  the  three  initial  face  map  parameters  (by 
calling  the  routine  ProjectPCA.m  described  below), 

•  facemapA.m  which  computes  the  range  of  possible  sub-images  contained 
in  the  image  provided  as  an  input,  and  leads  to  the  best  extracted  face  (by 
calling  the  routine  ProjectPCA.m  described  below),  FFT  implementation. 

•  facemap.m  which  duplicates  the  functions  facemapini.m  and  facemapA.m, 
without  fft  implementation  (much  slower) 

Manual  face  location  extraction  is  implemented  with  the  functions: 

•  cropimage.m  which  loads  the  requested  image  in  MATLAB,  calls  the 
function  cropface.m  which  returns  the  cropped  image,  where  the  cropped 
is  done  manually  by  the  operator. 

PCA-based  dimension  reduction  schemes  is  implemented  with  the  function: 

•  pea. m  which  returns  the  specified  number  of  projection  directions. 
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The  LDA-based  dimension  reduction  scheme  is  implemented  with  the  functions: 


•  pca.m  which  is  first  applied  to  address  the  singularity  issue  discussed  ear¬ 
lier  in  Chapter  IV, 

•  fld.m  which  returns  projection  directions  which  best  discriminate  the  data. 

Class-specific  centroid  computation  is  implemented  with  the  function: 

•  meanclass.m  which  takes  the  projection  direction  information  obtained 
with  pca.m  or  fld.m  and  computes  the  centroids  of  each  class  selected  by 
the  operator. 

Classification  perfonnance  is  computed  with  the  function: 

•  classify. m  which  uses  as  inputs  either  the  projection  direction  information 
or  the  images,  compares  to  the  class-specific  centroids,  and  returns  a  class 
decision. 

Cross-validation  results  are  implemented  with  the  two  functions: 

•  CrossValidateConstant.m  that  generates  the  error  rate  histograms  for  a 
60%  training/40%  testing  data  split,  where  the  number  of  images  per  class 
is  kept  constant.  (Results  provided  in  Chapter  V  use  this  routine). 

•  CrossValidate.m  that  generates  the  error  rate  histograms  for  60%  training 
40%  testing  data  split,  where  the  number  of  images  per  class  is  selected 
randomly  and  is  NOT  kept  constant. 


Working  routines: 

•  ProjectPCA.m  loads  a  set  of  user-specified  cropped  images,  computes  the 
PCA-based  projection  directions,  save  all  relevant  information  in  the  file 
face.dat.  This  routine  also  computes  the  error  rate  as  a  function  of  the 
number  of  eigenvectors,  as  shown  in  Figure  15,  Chapter  V. 

•  sortem.m  sorts  the  eigenvector  matrix  by  decreasing  eigenvalue  value. 
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•  InfraredTotalExitance.m  generates  the  total  black  body  exitance  for  user- 
specified  wavelength  and  temperature  ranges,  used  to  plot  Figure  1,  Chap¬ 
ter  III. 

•  distance.m  computes  the  distances  between  two  matrices  on  a  column-by¬ 
column  basis. 

•  projectLDA.m  loads  cropped  image  in  MATLAB,  applies  the  LDA  dimen¬ 
sion  reduction  scheme  for  a  user-specified  dataset  (training/testing  split  is 
user-specified),  investigates  further  dimension  reduction,  and  returns  the 
classification  performance.  Also  offers  the  option  to  load  the  original  im¬ 
ages. 
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%%%0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o%%% 


%  function  [W,m,  Amean, EVA]=pca(A,n) 

%  This  function  computes  the  pca(Principal  Component  Analysis)  of  the  data  contained 
%  in  the  array  A. 

%  input:  A  (k  x  n)  contains  the  data  in  each  column,  k  =  dimension  of  data  vector 
%  n  =  number  of  data  samples 

%  output:  W  contain  n  eigenvectors,  one  in  each  column 
%  m-  fraction  of  the  variance  in  n  eigenvectors 

%  Amean  -  mean  of  the  data  contained  in  A  in  a  column 

%  Ad  -  data  less  the  mean  Amean,  all  the  data  in  columns 
%  EVA  contain  the  eigenvalues  corresponding  to  the  eigenvectors  of  W  on  the 
%  diagonal 

%  Obs.  The  eigenvectors  associated  with  the  eigenvalues  smaller  that  1/1000  of  the 
%  largest  eigenvalue  are  eliminated 
%  Diogo  C.  Pereira  1LT  BRAF 
%  07/26/02 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


function  [W,m, Amean, Ad, EVA]=pca(A,n); 

[Adimk,Adimn]=size(A); 
if  (nargin==l) 
n=Adimn; 
end 

if  (n>Adimn) 
n=Adimn; 
end 

%  Amean  is  the  mean  of  A  using  the  columns  as  elements 
Amean=mean(  A, 2) ; 

%  Ad  is  the  difference  between  A  and  Amean 
Ad=A-Amean*ones(  1  ,Adimn); 

%eigenvectors  (columns  of  Vectors)  and  eigenvalues  (diag  of  Values) 

[EVE, EVA]  =  eig(Ad'*Ad); 

%  obtain  index  of  eigenvalues  greater  than  (0.001  times  greater  eigenvalue) 
[I]=find(EVA>(0.00 1  *max(max(EVA)))); 

EVAC=zeros(size(EVA)); 

EVAC(I)=EVA(I);  %EVAC  (Eigenvalue  Conditioned)  will  have  only  values  greater 
than  a  minimum  value  otherwise  the  value  is  set  to  zero 
%  Obtain  n  eigenvectors 

%EVAINV  is  the  matrix  containing  the  inverse  of  the  eigenvalues  on  the  diagonal,  it 
contains  zero  if  the  eigenvalue  was  zero 
%  and  the  inverse  of  the  eigenvalue  if  the  element  was  larger  than  zero 
EVAINV=zeros(Adimn);  %  EVAIN  has  size  n  x  n  since  n  is  the  maximum  number  of 
eigenvalues 

EVAINV(I)=E/(sqrt(EVA(I))); 

U=A*EVE*EVAINV ; 
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%Sort  the  vectors/values  according  with  the  absolute  value  of  the  eigenvalue  EVAC  and 
eliminates  the  columns  corresponding 
%  to  zero  eigenvalues 
[W,EVA]=sortem(U,  EVAC); 
if  (n>size(W,2)) 
n=size(W,2); 
end 

W=W(:,l:n); 
temp=diag(EVA); 
m=sum(temp(  1  :n))/sum(temp); 
return 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


%  function  [W,D]=fld(A,C) 

%  This  function  computes  the  Linear  Discriminant  W  that  maximizes 
%  the  ratio  |Wt  x  SB  x  W|/|Wt  x  SW  x  W|  where  SW  is  the  within  class 
%  scatter  matrix,  and  SB  is  the  between  class  scatter  matrix 
% 


%  input:  A-  data  matrix,  containing  data  in  columns 
%  C-row  array  containing  numbers  representing  the  classes  of 

%  the  elements  in  A 

%  output:  W-  weight  matrix  contains  the  vectors  in  columns 
%  D  -  matrix  containing  the  eigenvalues  on  the  diagonal 

%  last  modified:  08/19/2002 
%  Diogo  Pereira 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


function  [W,D]=fld(A,C); 

MeanG=mean(A,2);%MeanG  is  the  general  mean  of  the  data  in  A 
Ad=A-MeanG*ones(l,size(A,2));  %Ad  is  matrix  A  minus  the  mean  of  A 
ST=Ad*Ad';%  ST  is  the  total  scatter  matrix 
clear  Ad 

Nclass=0;  %Nclass  counts  the  number  of  classes  on  A 

%Compute  the  variance  within  all  the  classes  %%%%%%%%%%%%%%%%%%% 
SW=zeros(size(A,l));  %SW  is  the  within  class  scatter  matrix 
while  size(C,2)>0, 
w=C(l,l); 

[I]=fmd(C==w); 

MeanClass=mean(  A( :  ,I),2) ; 

Ad=A( :  ,1)-  MeanClass  *  ones(  1  ,size( A( :  ,I),2)); 

SW=SW+Ad*Ad’;  %adds  the  within  class  scatter  matrix  of  each  class 
A(:,I)=[];  %eliminates  the  data  already  used 

C(I)=[];  %eliminates  the  number  corresponding  to  the  class  from  C  vector 
Nclass=Nclass+l; 
end 
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%%%0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o%%% 


N=min(size(A,l),  Nclass-1);%N  will  be  the  minimum  between  the  dimension  of  the  data 
and  Nclass-1 

[W,EVA]=eigs(SB,SW,N); 

%Set  to  zero  the  eigenvalues  not  equal  to  a  finite  value 

[I]=find(isfinite(EVA)==l);  %Obtain  the  index  of  the  eigenvalues  that  are  finite 
EVAC=zeros(size(EVA)); 

E  V  AC  (I )=E  V  A(I) ; 

EVA=EVAC;%EVA  will  contain  just  the  eigenvalues  that  are  finite 
%  obtain  index  of  eigenvalues  greater  than  0.001  times  greater  eigenvalue 
[I]=find(abs(EVA)>(0.00 1  *max(max(abs(EVA))))); 

EVAC=zeros(size(EVA)); 

EVAC(I)=EVA(I);  %EVAC  will  have  only  values  greater  than  a  minimum  value  other¬ 
wise  the  value  is  set  to  zero 

[W,D]=sortem(W,EVAC);%  order  and  eliminates  eigenvector  corresponding  to  zero  ei¬ 
genvalue 

%  this  part  nonnalizes  W  so  that  the  norm  of  each  column  will  be  one 
NW=ones(size(W,  1 ),  1  )*sqrt(sum((W.  A2),  1 )); 

W=W./NW; 

return 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


%  CrossValidateConstant 

%  This  programs  test  the  2  schemes  using  cross  validation  6  training/4  testing 
%  The  number  of  elements  in  each  class  is  left  constant 
%  Diogo  Pereira 


%  12/14/2002 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


clear  all; 
close  all; 

%  Load  the  images  %%%%%%%%%%%%%%%%%%%%%% 

A=[];  %  A  contains  images  in  columns 

T=[];  %  T  contains  the  class  number  of  each  picture 

Crop- a'; 

N_pictures=10;  %  number  of  pictures  for  each  person  number 

epct=0;  %  accumulator  of  the  number  of  errors  per  class 

epct  100=0; 

epct40=0; 

epct40Wl=0; 

epct40W2=0; 

epct40W3=0; 

epct40W4=0; 

epct40W5=0; 

epctLDA=0; 

Person=[l  2345689  11  12  13  14  15  16];%  Person  numbers  contained  on  the  files 
Section=l; 
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[A,T]=LoadImage(A,T, Person,  N_pictures, Section, Crop); 

Section=5; 

[A, T]=LoadImage(A,T, Person,  N_pictures, Section, Crop); 

Section=6; 

[A, T]=LoadImage(A,T, Person,  N_pictures, Section, Crop); 

%  Obtain  the  samples  to  be  used  on  the  testing  and  remove  it  from  the  training  images 

%  A  will  contain  the  training  images 

%  B  will  contain  the  testing  images 

%  T  is  the  class  number  of  the  images  on  A 

%  TB  is  the  class  number  of  the  images  on  B 

Atemp=A; 

Ttemp=T; 

N_trials=300; 

N_samples=4*14*3;  %  we  are  using  6/4  cross  validation 

for  j=l  :N_trials 
A=Atemp; 

T=Ttemp; 

N_Person=14; 

B=[]; 

TB=[]; 

Remove=[]; 

%  Generates  testing  and  training  sets 

for  i=l:(N_Person*3)  %  i  goes  from  1  to  the  total  number  of  pictures 
f=randperm(10); 

Ll=(i-l)*10+l; 

L2=i*10; 

A(:,Ll:L2)=A(:,(i-l)*10+f);  %  mix  the  samples  in  each  class 
T(:,Ll:L2)=T(:,(i-l)*10+f);  %  adjust  the  class  numbers 
IndSamples=[Ll  Ll+1  L1+2L1+3]; 

B=[B  A(:,IndSamples)]; 

TB=[TB  T(:,IndSamples)]; 

Remo ve= [Remove  IndSamples]; 
end 

A(:,Remove)=[]; 

T(:,Remove)=[]; 


%  Compute  the  pea  of  the  training  images  A 
[W,m,Amean,Ad]=pca(A); 

%  Compute  the  projection  matrix  P 
P=W'*(Ad); 

%%%%  used  to  do  cross  validation  directly  without  pea  or  Ida 
PD=A; 
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%  AmeanD=mean(  A, 2) ; 

PD=A-Amean*ones(  1  ,size(A,2)); 

%  Compute  the  centroid  of  each  class  contained  in  Person 
C=meanclass(P,T,  Person); 

CD=meanclass(PD,T,  Person); 

%C(1, :)=[]; 

%  Using  the  testing  images  B 
Ad=B-Amean*ones(  1  ,size(B,2)); 

%  Compute  the  projection  matrix  P 
P=W'*(Ad); 

%P(1, :)=[]; 

%  Classify  the  testing  images 
dmin=100000000; 

%  classifying  the  images  directly 
PD=Ad; 

[D]=classify(PD,CD,dmin); 

NZD=fmd(D~=0); 

D(NZD)=Person(D(NZD)); 

error=TB-D; 

n_erxorDirect(j)=sum((error~=0),2); 

%computes  the  number  of  errors  per  class 

errornz=(error~=0); 

epc=reshape(erromz,4,42); 

epcpl=epc(:,l:14); 

epcp2=epc(:,  1 5 :28); 

epcp3=epc(:,29:42); 

epc=[epcp  1  ;epcp2;epcp3] ; 

epc=sum(epc,l); 

epct=epct+epc; 

%  limiting  the  size  to  Neig=number  of  eigenvectors 

%  using  100  eigenvectors 

Neig=100; 

N  eig=min(N  eig,  length(C)) ; 

Templ=P(l:Neig,:); 

Temp2=C(l:Neig,:); 

[D]=classify(Templ,Temp2,dmin); 

NZD=fmd(D~=0); 

D(NZD)=Person(D(NZD)); 

error=TB-D; 

n_errorPCA  1 00(j)=sum((error~=0),2); 

%computes  the  number  of  errors  per  class 

errornz=(error~=0); 

epc=reshape(erromz,4,42); 

epcpl=epc(:,l:14); 

epcp2=epc(:,  15:28); 

epcp3=epc(:,29:42); 
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epc=[epcp  1  ;epcp2;epcp3] ; 

epc=sum(epc,l); 

epct  1 00=epct  1  OO+epc; 

%using  40  eigenvectors 
Neig=40; 

Templ=P(l:Neig,:); 

Temp2=C(l:Neig,:); 

[D]=classify(Templ,Temp2,dmin); 

NZD=find(D~=0); 

D(NZD)=Person(D(NZD)); 

%TB 

%D40=D 

error=TB-D; 

n_errorPCA40(j)=sum((error~=0),2); 

%computes  the  number  of  errors  per  class 

erromz=(error~=0); 

epc=reshape(erromz,4,42); 

epcpl=epc(:,l:14); 

epcp2=epc(:,  15:28); 

epcp3=epc(:,29:42); 

epc=[epcp  1  ;epcp2;epcp3] ; 

epc=sum(epc,l); 

epct40=epct40+epc; 

%using  40  eigenvectors  without  the  first  one 
Neig=40; 

Templ=P(2:Neig+l,:); 

Temp2=C(2:Neig+ 1 , 

[D]=classify(Templ,Temp2,dmin); 

NZD=fmd(D~=0); 

D(NZD)=Person(D(NZD)); 

error=TB-D; 

n_errorPCA40W  1  (j)=sum((error~=0),2); 

%computes  the  number  of  errors  per  class 

errornz=(error~=0); 

epc=reshape(erromz,4,42); 

epcpl=epc(:,l:14); 

epcp2=epc(:,  15:28); 

epcp3=epc(:,29:42); 

epc=[epcp  1  ;epcp2;epcp3] ; 

epc=sum(epc,l); 

epct40Wl=epct40Wl+epc; 

%using  40  eigenvectors  without  the  first  two  eigenvectors 
Neig=40; 

Templ=P(3:Neig+2,:); 

Temp2=C(3:Neig+2,:); 

[D]=classify(Templ,Temp2,dmin); 
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NZD=find(D~=0); 

D(NZD)=Person(D(NZD)); 

error=TB-D; 

n_errorPCA40W2(j)=sum((crror~=0),2); 

%computes  the  number  of  errors  per  class 

erromz=(error~=0); 

epc=reshape(erromz,4,42); 

epcpl=epc(:,l:14); 

epcp2=epc(:,  1 5 :28); 

epcp3=epc(:,29:42); 

epc=[epcp  1  ;epcp2;epcp3] ; 

epc=sum(epc,l); 

epct40W2=epct40W2+epc; 

%using  40  eigenvectors  without  the  first  three  eigenvectors 
Neig=40; 

Templ=P(4:Neig+3,:); 

Temp2=C(4:Neig+3 , 

[D]=classify(Temp  1  ,Temp2,dmin); 

NZD=find(D~=0); 

D(NZD)=Person(D(NZD)); 

error=TB-D; 

n_errorPCA40W3(j)=sum((error~=0),2); 

%computes  the  number  of  errors  per  class 

errornz=(error~=0); 

epc=reshape(erromz,4,42); 

epcpl=epc(:,l:14); 

epcp2=epc(:,  15:28); 

epcp3=epc(:, 29:42); 

epc=[epcp  1  ;epcp2;epcp3] ; 

epc=sum(epc,l); 

epct40W3=epct40W3+epc; 

%using  40  eigenvectors  without  the  first  four  eigenvectors 
Neig=40; 

Templ=P(5:Neig+4,:); 

Temp2=C(5:Neig+4,:); 

[D]=classify(Templ,Temp2,dmin); 

NZD=find(D~=0); 

D(NZD)=Person(D(NZD)); 

error=TB-D; 

n_errorPCA40W4(j)=sum((error~=0),2); 

%computes  the  number  of  errors  per  class 

errornz=(error~=0); 

epc=reshape(erromz,4,42); 

epcpl=epc(:,l:14); 

epcp2=epc(:,  15:28); 

epcp3=epc(:,29:42); 
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epc=[epcp  1  ;epcp2;epcp3] ; 

epc=sum(epc,l); 

epct40W4=epct40W4+epc; 

%using  40  eigenvectors  without  the  first  five  eigenvectors 
Neig=40; 

Templ=P(6:Neig+5,:); 

Temp2=C(6:Neig+5,:); 

[D]=classify(Templ,Temp2,dmin); 

NZD=fmd(D~=0); 

D(NZD)=Person(D(NZD)); 

error=TB-D; 

%D40W5=D 

n_errorPCA40W5(j)=sum((error~=0),2); 

%computes  the  number  of  errors  per  class 

erromz=(error~=0); 

epc=reshape(erromz,4,42); 

epcpl=epc(:,l:14); 

epcp2=epc(:,  1 5:28); 

epcp3=epc(:,29:42); 

epc=[epcpl;epcp2;epcp3]; 

epc=sum(epc,l); 

epct40W5=epct40W5+epc; 

0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o%%% 

%  LDA 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


%  Compute  the  pea  of  the  training  images 
temp=(N_pictures-l)*length(Person); 

[W,m,Arnean,Ad,EVA]=pca(A,temp);  %  This  reduces  the  dimension  to  N-c  =  36  -  6=30 
the  best  result  was  with  29 
%  Compute  the  projection  matrix  P 
P=W'*(Ad); 

%  P(l, :)=[]; 

%  Compute  the  projection  matrix  P  obtained  with  the  FLD(Fisher  Finear  Discriminant) 
[Wopt,D]=fld(P,T); 

P=Wopt' *P; 

%  Computing  the  centroid  of  each  class 
C=meanclass(P,T,Person); 

%  Using  the  testing  data 
Ad=B-Amean*ones(  1  ,size(B,2)); 

P=W’*(Ad); 

%P(1, :)=[]; 

P=Wopf*P; 

%  Classify  the  testing  images 
dmin=100000000; 

[D]=classify(P,C,dnhn); 

NZD=find(D~=0); 
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D(NZD)=Person(D(NZD)); 

error=TB-D; 

n_errorLD  A(j  )=sum((error~=0 )  ,2) ; 

%computes  the  number  of  errors  per  class 

erromz=(error~=0); 

epc=reshape(erromz,4,42); 

epcpl=epc(:,l:14); 

epcp2=epc(:,  15:28); 

epcp3=epc(:, 29:42); 

epc=[epcp  1  ;epcp2;epcp3] ; 

epc=sum(epc,l); 

epctLDA=epctLDA+epc; 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


end 

h=figure; 

L 1  =min(n_errorDirect); 

L2=max(n_errorDirect); 
temp  1 =linspace(0,L2,L2+ 1 ); 

N 1  =hist(n_errorDirect,temp  1 ); 

plot(templ/N_samples*100,Nl/N_trials*100,’bd-'); 
title('Histogram  of  the  number  of  errors'); 
xlabel('Percentage  of  errors'); 
ylabel('Percentage  of  the  simulations'); 
grid; 

saveas(h,'CrossValidationDirect6_4.fig’); 

h=figure; 

L 1  =min(n_errorPCA  100); 

L2=max(n_errorPC  A 1 00); 
temp  1 =bnspace(0,L2,L2+ 1 ); 

N 1  =hist(n_errorPCA  1 00,temp  1 ); 
plot(templ/N_samples*100,Nl/N_trials*100,’bd-'); 
title('Histogram  of  the  number  of  errors'); 
xlabel('Percentage  of  errors'); 
ylabel('Percentage  of  the  simulations'); 
grid; 

saveas(h,'CrossValidationl00eig_6_4.fig'); 

h=figure; 

L 1  =min(n_errorPC  A40); 

L2=max(n_errorPCA40); 
temp  1 =bnspace(0,L2,L2+ 1 ); 

N  l=hist(n_errorPCA40,temp  1); 
plot(templ/N_samples*100,Nl/N_trials*100,’bd-'); 
title('Histogram  of  the  number  of  errors'); 
xlabel('Percentage  of  errors'); 
ylabel('Percentage  of  the  simulations'); 
grid; 
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saveas(h,'CrossValidation40eig_6_4.fig'); 

h=figure; 

L  l=min(n_errorPCA40Wl); 
L2=max(n_errorPCA40W  1 ); 
temp  1 =linspace(0,L2,L2+ 1 ); 

N 1  =hist(n_errorPCA40W  1  ,temp  1 ); 
plot(templ/N_samples*100d41/N_trials*100,’bd-'); 
title('Histogram  of  the  number  of  errors'); 
xlabel('Percentage  of  errors'); 
ylabel('Percentage  of  the  simulations'); 
grid; 

saveas(h,'CrossValidation40Wleig_6_4.fig'); 

h=figure; 

L  l=min(n_errorPCA40W2); 
L2=max(n_errorPCA40W2); 
temp  1 =linspace(0,L2,L2+ 1 ); 

N  l=hist(n_errorPCA40W2,temp  1); 
plot(templ/N_samples*100,Nl/N_trials*100,’bd-'); 
title('Histogram  of  the  number  of  errors'); 
xlabel('Percentage  of  errors'); 
ylabel('Percentage  of  the  simulations'); 
grid; 

saveas(h,'CrossValidation40W2eig_6_4.fig'); 

h=figure; 

L  l=min(n_errorPCA40W3); 
L2=max(n_errorPCA40W3); 
temp  1 =linspace(0,L2,L2+ 1 ); 

N 1  =hist(n_errorPCA40W  3  ,temp  1 ); 
plot(templ/N_samples*100,Nl/N_trials*100,’bd-'); 
title('Histogram  of  the  number  of  errors'); 
xlabel('Percentage  of  errors'); 
ylabel('Percentage  of  the  simulations'); 
grid; 

saveas(h,'CrossV  alidation40  W3  eig_6_4 .  fig'); 
h=figure; 

L  l=min(n_errorPCA40W4); 
L2=max(n_errorPCA40W4); 
temp  1 =linspace(0,L2,L2+ 1 ); 

N  l=hist(n_errorPCA40W4,temp  1); 
plot(templ/N_samples*100,Nl/N_trials*100,’bd-'); 
title('Histogram  of  the  number  of  errors'); 
xlabel('Percentage  of  errors'); 
ylabel('Percentage  of  the  simulations'); 
grid; 

saveas(h,'CrossValidation40W4eig_6_4.fig'); 
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h=figure; 

L  l=min(n_errorPCA40W5); 

L2=max(n_errorPCA40W5); 
temp  1 =linspace(0,L2,L2+ 1 ); 

N 1 =hist(n_errorPCA40W  5  ,temp  1 ); 
plot(templ/N_samples*100,Nl/N_trials*100,’bd-'); 
title('Histogram  of  the  number  of  errors'); 
xlabel('Percentage  of  errors'); 
ylabel('Percentage  of  the  simulations'); 
grid; 

saveas(h,'CrossV  alidation40  W5  eig_6_4 .  fig'); 

%  hold  on 
h=figure; 

L 1  =min(n_errorLD  A); 

L2=max(n_errorLD  A) ; 
temp2=linspace(0,L2,L2+l); 

N  2=hist(n_errorLD  A,  temp2 ) ; 
plot(temp2/N_samples*  100,N2/N_trials*  100,’rd-'); 
title('Histogram  of  the  number  of  errors'); 
xlabel('Percentage  of  errors'); 
ylabel('Percentage  of  the  simulations'); 
grid; 

saveas(h,'CrossV  alidationLD  A_6_4 .  fig'); 

%  saving  the  error  per  class 

ns=((N_samples/ 1 4)*N_trials)/ 100; 

epct=epct/ns; 

epct  1 00=epct  1 00/ns; 

epct40=epct40/ns ; 

epct40Wl=epct40Wl/ns; 

epct40W2=epct40W2/ns; 

epct40W3=epct40W3/ns; 

epct40W4=epct40W4/ns; 

epct40W5=epct40W5/ns; 

epctLDA=epctLDA/ns; 

M=[Person;epct;epctl00;epct40;epct40Wl;epct40W2;epct40W3;epct40W4;epct40W5;ep 

ctLDA]; 

filename- ErrorPerClass6_4perc'; 

wk  1  write(filename,M,4,2) ; 

ns=  1 00/((N_samples/ 1 4)  *N_trials); 

epct=epct/ns; 

epct  1 00=epct  1 00/ns; 

epct40=epct40/ns ; 

epct40Wl=epct40Wl/ns; 

epct40W2=epct40W2/ns; 

epct40W3=epct40W3/ns; 

epct40W4=epct40W4/ns; 
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epct40W5=epct40W5/ns; 

epctLDA=epctLDA/ns; 

M=[Person;epct;epctl00;epct40;epct40Wl;epct40W2;epct40W3;epct40W4;epct40W5;ep 

ctLDA]; 

filename='ErrorPerClass6_4abs'; 
wk  1  write(filename,M,4,2) ; 


%  LocateClassify 

%  this  routine  locates  a  face  and  classify  it  according  with  the  PCA  method 

% 


%  input:  name  of  the  file  containing  the  image 

%  W  eigenvectors,  H  high  of  the  image  on  the  database,  C  centroids  of  the  classes 

%  Amean  mean  of  the  data  used  to  obtain  the  eigenvectors 

%  output:  number  of  the  class  the  object  was  classified  to 

%  Obs.:  You  need  to  run  the  program  ProjectPCA  and  ProjectLDA  first  to  obtain 

%  the  files  faces  and  faces2 

%  Diogo  Pereira 

%  last  update:  10/09/2002 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


close  all 
clc 

%  obtaining  the  name  of  the  image 
name=input('Type  the  name  of  the  file  ->','s’); 

T=3; 

B=readpgm8(name);  %  read  pgm  files 
[HI,WI]=size(B);%  obtaining  the  size  of  the  image 

%  displays  the  image%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

figure('Units','pixels','Position',[100  100  WI  HI]); 
imagesc(B) 

set(gca, 'Units', 'pixels', 'Position', [0  0  WI  HI]); 

colonnap(gray) 

pause(.8) 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


disp( ['Locating  the  face,  wait...’]) 

%  Locates  and  obtains  the  face%%%%%%%%%%%%%%%%%%%%%%%%%% 

B=locatehead(B);  %locates  the  head 

B=double(B); 

load  faces  W  H  Amean  C 

load  faces2  W2  C2 

[m,n]=size(W); 

WI=m/H; 

F=facemap(W,H, Amean, B);  %  returns  F  the  facemap 
[M,I]=min(F); 

[N,yini]=min(M); 
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xini=I(yini); 

face=B(xini:(xini+H-l),yini:(yini+WI-l));  %  obtain  the  face  on  the  image 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


%  displays  the  face  image%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

figure('Units', 'pixels', 'Position',  [400  100  WI  H]); 
imagesc(face) 

set(gca, 'Units', 'pixels', 'Position', [0  0  WI  H]); 
colonnap(gray) 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


%  classify  the  face  image%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

face=reshape(face,WI*H,  1); 

face=face-Amean; 

%classify  using  the  PCA  scheme 

PB=W'*face; 

dmin=100000000; 

kmax=size(PB,l); 

k=kmax; 

D=classify(PB(l  :k,:),C(  1  :k,:),dmin); 

disp(['You  belongs  to  class  ’,num2str(D),'  using  the  PCA  classification  scheme’]) 

%classify  using  the  LDA  scheme 

PC=(W2)'*face; 

dmin=100000000; 

kmax=size(PC,l); 

k=kmax; 

D=classify(PC(l  :k,:),C2(  1  :k,:),dmin); 

disp(['You  belongs  to  class  ’,num2str(D),'  using  the  LDA  classification  scheme’]) 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

%  Function  [A, T]=LoadImage(A,T, Person,  N_pictures, Section); 

%  This  function  loads  images  to  the  matrix  A 

%  input:  A  -  matrix  containing  images  on  each  column,  may  be  empty 
%  T  -  row  containing  the  class  number  of  each  picture,  may  be  empty 
%  N_pictures  -  number  of  pictures  in  each  person  number 
%  Section  -  photo  section  number 
%  Crop-  letter  indicating  the  crop  section 
%  Diogo  Pereira 
%  12/14/2002 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


function  [B,T]=LoadImage(  A, T, Person,  N_pictures, Section, Crop); 
N_class=length(Person); 
for  i=l:N_class 
for  j= 1  :N_pictures 

imnum  I  =  num2str(Person(i)); 
im_num2=  num2str(j); 
s=num2str(  S  ection) ; 

img_name  =  strcat(im_numl,'-’,im_num2,'-',s,'-’,Crop); 


72 


img=  imread(img  name, ’bmp'); 

[diml,dim2]  =  size(img); 
x=reshape(img,diml  *dim2, 1); 

A=[A  x];  %  A  contains  one  image  in  each  row 
T=[T  Person(i)];%  T  contains  the  class  of  each  image 
end 
end 


clear  x; 

B=double(A); 

Return 


%  function  [C]=meanclass(A,T) 

%  this  function  computes  the  mean  of  the  coordinates  of  pictures  of  each  class 
%  input :  A  -  matrix  containing  the  coordinates  of  each  picture 
%  in  columns 


%  T  -  vector  containing  the  class  number  of  each  picture 

%  P  -  vector  containing  the  class  numbers  that  we  want  to  obtain  the  mean 

%  output:  C  -  matrix  containing  the  mean  of  the  coordinates  of  the 
%  pictures  of  each  class  contained  in  T 

%  Diogo  Pereira  11/13/02 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


function  [C]=meanclass(A,T,P) 

C=[]; 

for  i=l  :length(P) 

Temp=P(i); 

TMean=find(T==Temp); 
Mean=mean(  A( :  ,TMean),2) ; 
C=[C  Mean]; 
end 


return 


%  Function  FacemapA  computes  the  square  of  the  difference 
%  between  a  subimage  and  its  projection  on  the  facespace; 

%  the  subimage  has  the  size  of  the  eigenface 

%  input:  W-  array  (m  x  n)  containing  the  eigenfaces  in  columns 

%  m  dimension  of  each  eigenface,  n  number  of  eigenfaces 

%  diml  -  dimension  1  of  each  eigenface,  dim2=m/diml 

%  Amean-  vector  containing  the  mean  of  the  data  used  to  create  the  eigenfaces 

%  A-  array  (p  x  q)  containing  the  image  to  be  analyzed 

%  output  C-  array  (r  x  t)  containing  the  difference  between  the  subimage  and  its  projec¬ 
tion  squared 

% 

%  Diogo  Pereira 
%  10/15/02 
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%%0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o%%% 


function  [C]=facemapA(W,dhnl  ,Amean,A,T  1  ,T2,T3) 
[M,N]=size(A); 

[P,L]=size(W); 

AF=FFT2(A,256,256);  %  AF  is  the  FFT  of  the  image 
Ll=M-diml+l; 

L2=N-(P/dim  1 )+ 1 ; 

T5=zeros(Ll,L2); 
for  i=l:L 

U=reshape(W  ( :  ,i),dim  1  ,P/< dim  1 ) ; 

UF=AF.*conj(FFT2(U,256,256)); 

U=real(IFFT2(UF)); 

U=U(1:L1,1:L2); 

U=U-T3(l,i); 

U2=U.*U; 

T5=T5+U2; 

end 

Amean=reshape(Amean,dim  1 ,  P/dim  1 ); 

AmeanF=FFT2(Amean,256,256); 

T4F=AF.*conj(AmeanF); 

T4F =real(IFFT2  (T  4F)) ; 

T4=2*T4F(1:L1,1  :L2); 

C=T1+T2-T4-T5; 

return 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


%  function  [Tl,T2,T3]=facemapini(A,B,W) 

%  This  function  computes  tenns  Tl,  T2  and  T3  of  the  five  terms  used  on  the  face  map 
%  equation 

%  input:  A-  (M  x  N)Image  where  the  face  is  located 
%  B-  Mean  face  of  the  eigenfaces  in  a  column 
%  W-  (K  x  L)Matrix  containing  eigenfaces  in  columns 
%  diml-  vertical  dimension  of  the  eigenface 
%  Diogo  Pereira 
%  last  update:  10/15/02 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


function  [T 1  ,T2  ,T3  ]=facemapini( A,B , W,dim  1 ); 
[M,N]=size(A); 

U=diml; 

P=length(B)/diml ; 

[K,L]=size(W); 

Ll=M-diml+l; 

L2=N-P+1; 

Tl=zeros(Ll,L2); 

A2=A.*A; 

%  Computing  T 1  using  FFT2 
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A2F=FFT2(A2,256,256); 
WindowF=conj'(FFT2(ones(diml,P),256,256)); 
T 1  =IFFT2(A2F.  *  WindowF); 
Tl=real(Tl(l:Ll,l:L2)); 

%  Computing  T2 
B=reshape(B,  1 ,11*?); 

T2=B*B’; 

T2=ones(Ll,L2)*T2; 

%  Computing  T3 

T3=(B*W); 

return 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


%  function  [B]=locatehead(A) 

%  this  function  will  locate  the  boundaries  of  the  head 
%  and  will  crop  the  image  around  the  head 
%  it  locates  just  the  horizontal  boundaries 
%  input:  A  -  original  image 
%  output:  B  -  cropped  image 
% 


%  Diogo  Pereira 
%  last  update:  10/1 1/02 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


function  [B]=locatehead(A); 

A=double(A); 

[M,N]=size(A); 

P=l; 

F=N; 

Z=zeros(M,l); 

D=[A  Z]-[Z  A]; 
d=floor(M/6); 

Pl=fmd(D(2*d,:)>25); 

P2=find(D(3*d,:)>25); 

P3=find(D(4*d,:)>25); 

%  computing  P  the  minimum  among  P 1 ,  P2  and  P3 
if  ~iscmpty(P  1 ) 

p=pi(1); 

end 

if  ~iscmpty(P2) 

P=min([P  P2(  1 )]); 
end 

if  ~isempty(P3) 

P=min([P  P3(  1 )]); 
end 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


Fl=find(D(d*2,:)<-20); 
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F2=find(D(d*3,:)<-20); 

F3=find(D(d*4,:)<-20); 

%  computing  F  the  maximum  among  FI,  F2  and  F3 
if  -isempty(Fl) 

F=Fl(length(Fl)); 

end 

if  ~isempty(F2) 

F=max([F  F2(length(F2))]); 
end 

if  ~isempty(F3) 

F=max([F  F3(length(F3))]); 
end 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


B=A(:,P:F); 

Return 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

%  function  [D]=classify(PB,P,dmin) 

%  This  function  classify  PB  in  classes,  each  column  of  P  represents  a  class 
%  dmin  is  the  minimum  distance  from  each  class.  If  the  distance  of  PB  to  the 
%  corresponding  class 

%  is  larger  than  dmin,  the  anwswer  return  class  O.  The  distance  is  the  Euclidian  distance 
%  input:  PB  matrix  containig  the  data  to  be  classified,  each  column  contain  the  data 
%  P  matrix  containing  the  coefficient  of  each  class,  each  column  contain  the  data  of 

%  one  class 

%  dmin  minimum  distance  from  each  class,  it  is  a  scalar 

%  output:  D  row  vector  containing  the  number  of  the  class  that  each  object  was  classified 
%  to 
% 


%  Diogo  Pereira  07/28/2002 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


function  [D]=classify(PB,P,dmin); 
[L,M]=size(PB);  %there  are  M  objects 
[N,K]=size(P);  %there  are  K  classes 
D=zeros(l,M); 
for  i=l:M 


dif=P-PB(:,i)*ones(  1  ,K); 
[dist,pos]=min((  sum((dif).A2,l) )); 
if  dist<dmin 
D(l,i)=pos; 
end 
end 
return 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


%  function  [C]=distance(A,B) 
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%  this  function  computes  the  distance  between  the  elements  in  A  and  the  elements  in  B 

% 

%  input  A=  array  containing  vectors  in  columns 
%  B=  array  containing  vectors  in  columns 
%  output  C  =  array  containing  the  distance  between  the 
%  Diogo  Pereira 
%  11/01/02 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


function  [C]=distance(A,B); 
[L,M]=size(A);  %there  are  M  objects 
[N,K]=size(B);  %there  are  K  objects 
C=zeros(M,K); 
for  i=l:M 


dif=B  -  A( : ,  i)  *  ones(  1  ,K) ; 
C(i,:)=sqrt((  sum((dif).A2,l) )); 
end 
return 


function  [NV,ND]  =  sortem(V,D) 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


%[V,D]  =  SORTEM(V,D) 

%Sorts  the  columns  of  V  along  with  the  absolute  value  of  the  elements  of  D  and 
%  eliminates  the  column  of  V 
%  corresponding  to  zero  on  the  diagonal  of  D 


% 


%  Diogo  Pereira,  2002 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


dvec  =  diag(D);  %obtain  the  values  of  the  diagonal  and  insert  in  a  vector 
%eliminating  the  column  on  V  corresponding  to  zero  on  D 
[l]=find(dvec==0); 
dvec(I)=[]; 

V(:,I)=[]; 

NV  =  zeros(size(V)); 

%sort  the  elements  of  dvec  in  descending  order  according  with  the  absolute  value  of  dvec 
[L,index_dv]  =  sort(abs(dvec)); 
indexdv  =  fhpud(indexdv); 
dvec=dvec(index_dv); 

%insert  the  elements  of  dvec  on  the  diagonal 
ND=diag(dvec); 

%sort  the  columns  of  V  according  with  index  dv 
NV=V  ( :  ,index_dv); 

Return 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


%  ConverPGMtoBMP  -this  file  loads  all  the  images  stored  in  pgm  format  and  save 
%  the  cropped  image  in  bmp  format 
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% 

%  Diogo  Pereira 
%  12/14/2002 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


Person=[l  2345689  11  12  13  14  15  16];%  Person  numbers  contained  on  the  files 

N_pictures=10;  %  number  of  pictures  for  each  person  number 

Section=l; 

load  faces  W  H  Amean  C  %Amg 
for  i=l  dength(Person) 
for  j= 1  :N_pictures 

im  num  I  =  num2str(Person(i));  im_num2=  num2str(j);im_num3=num2str(Section); 
imgname  =  strcat(im_numl,'-',im_num2,'-',im_num3); 

A=readpgm8(img_name);  %  read  pgm  files 
%Head=locatehead(A);  %locates  the  head 
B=double(A); 

[m,n]=size(W); 

WI=m/H; 

[Tl,T2,T3]=facemapini(B,Amean,W,H); 

Q=facemapA(W,H,Amean,B,T  1  ,T2,T3); 

[M,I]=min(Q); 

[N,yini]=min(M); 

xini=I(yini); 

face=A(xini:(xini+H-l),yini:(yini+WI-l));  %  obtain  the  face  on  the  image 
namefinal=[img_name  ’-b.bmp']; 
imwrite(face,namefinal,’bmp');  %  save  the  files  as  bmp 
end 
end 


%  cropface  -returns  a  portion  on  the  image  selected 
%  dragging  a  windows  with  a  mouse,  the  size  of  the  window  is 
%  given  by  H  times  W 

% 

%  function  [B]=cropface(A,W,  H) 

%  input  A  =  original  image,  matrix  containing  the  image 
%  W  =  number  of  pixels  wide 
%  H  =  number  of  pixels  high 
%output  B  =  image  selected 
% 


%  Diogo  Pereira 
%  last  update:  10/03/2002 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


function  [B]=cropface(A,W,H); 

[diml,dim2]=size(A); 

figure('Units',’pixels','Position',[100  100  dim2  diml]); 
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imagesc(A); 

colormap(gray); 

set(gca, 'Units', 'pixels', 'Position', [0  0  dim2  diml]); 
waitforbuttonpress; 

point  1  =  get(gcf,’CurrentPoint’);  %  button  down  detected 
rect  =  [pointl(l,l)  pointl(l,2)  W  H]; 

[P]  =  dragrect(rect); 

L=P(1,1);  %L  is  the  Left 
B=P(1,2);  %B  is  the  Bottow 
xini=diml-B-H; 
yini=L; 

B= A(xini :  (xini+H- 1  ),y  ini :  (y  ini+ W  - 1 )) ; 

Return 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

%  cropimage  -  this  program  crops  the  images 
%  The  user  click  on  the  figure  and  a  window  is  dragged 
%  when  the  user  releases  the  mouse,  the  portion  contained 
%  within  the  rectangle  is  cropped  and  saved  with  the  same  name 
%  as  the  original  file  but  with  a  ’-a’  added  to  the  file  name 
%  input:  file  name,  typed 
%  output:  image  cropped 
% 


%  Diogo  Pereira 
%  last  update:  10/03/2002 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


clear 
close  all 
clc 


name-’; 

%  data  input 
while  name~='fim' 

name=input('Type  the  name  of  the  file  ->','s'); 
temp=name(  1,1:3); 
if  temp=='fim' 
break 
end 

close  all 

a=readpgm8(name);  %  read  pgm  files 
W=45; 

H=60; 

b=cropface(a,W,H); 

figure('Units',’pixels','Position',[100  100  W  H]); 
imagesc(b); 

set(gca, 'Units', 'pixels', 'Position', [0  0  W  H]); 
colormap(gray); 
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namefinal=[name  '-a']; 
imwrite(b  ,namefmal,  'bmp') ; 
name=name(  1,1:3); 
end 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


%  ProjectPCA 

%  Diogo  Pereira  07/28/02 

%  This  program  saves  the  file  faces  used  in  other  programs 

% 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


clear 
close  all 
clc 

%  Testing  the  routine  using  the  first  N_object  faces  of  each  class  for  training  and  the  last 

%  10-N_object  faces  for  testing 

N_class=16;  %N_Class  is  the  number  of  classes 

N_object=10;  %N_object  is  the  number  of  objects  in  each  class 

A=[]; 

T=[]; 

%  W  is  an  array  containing  the  eigenvectors  of  the  covariance  matrix  in  each  column 
%  Load  the  training  images  into  A  matrix 
Personal  2  3  4  5  6  8  9  1 1  12  13  14  15  16]; 

N_class=length(Person); 

N_object=10; 
for  i=l:N_class 
for  j=l:N_object 

imnum  1  =  num2str(Person(i)); 
im_num2=  num2str(j); 

imgname  =  strcat(im_numl,'-',im_num2,'-l-a'); 
img=  irnread(img_name,’bmp'); 

[diml,dim2]  =  size(img); 
x=reshape(img,diml  *dim2, 1); 

A=[A  x];  %  A  contains  one  image  in  each  row 
T=[T  Person(i)]; 
end 
end 

clear  x; 

N_pictures=N_class*N_object;  %N_pictures  is  the  total  number  of  pictures  on  the  A  ar¬ 
ray 

A=double(A); 

%show  up  to  100  images  contained  in  the  array  A 
figure/ 1) 

for  i=l:min([N_pictures  100]); 
subplot/ 10, 10, i) 

trial_face=reshape(A(:,i),diml,dim2); 
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imagesc(trial_face) 
colormap(gray); 
axis  off; 
end 

H=diml ; 

%  Compute  the  pea  of  the  data  matrix 
[W,m,Amean,Ad]=pca(A); 

%  Compute  the  matrix  P 
P=W'*(Ad); 

%  Computing  the  centroid  of  each  class=the  average  of  the  coefficients 

C=[]; 

for  i=l  :N_class 

U 1  =(i- 1 )  *N_obj  ect+ 1 ; 

U2=U  1  +N_obj  ect- 1 ; 

C(:,i)=mean(P(:,(U  1  :U2)),2); 
end 

save  faces  W  Amean  H  C  %  stores  some  data  on  the  file  faces 
%  Loading  the  testing  data 
%  Load  the  test  images  into  B  matrix 
B=[]; 

T=[]; 

N_class=14; 

Person=[l  2  3  4  5  6  8  9  1 1  12  13  14  15  16]; 

N_class=length(Person); 

N_object=10; 
for  i=l  :N_class 
for  j= 1:10 

imnum  1  =  num2str(Person(i)); 
im_num2=  num2str(j); 

imgname  =  strcat(im_numl,'-',im_num2,'-5-a'); 
img=  imread(img_name,’bmp'); 

[diml,dim2]  =  size(img); 
x=reshape(img,diml*dim2,l); 

B=[B  x];  %  B  contains  one  image  in  each  row 
T=[T  Person(i)]; 
end 
end 

clear  x; 

B=double(B); 

LB=size(B,2); 

B=B -Amean*  ones(  1  ,LB); 

PB=W'*B; 

%  Computing  the  classification 
dmin=100000000; 
kmax=size(PB,  1); 
for  k=l:kmax 
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[D]=classify(PB(l:k,:),C(l:k,:),dmin); 

error=T-D; 

n_error(k)=sum((error~=0)  ,2); 
end 

figure(2) 

plot(n_error,'+b') 

xlabel('Number  of  eigenvectors'); 

ylabel(’N umber  of  Errors'); 

title('Classification  errors  using  140  test  faces'); 

grid; 

end 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


%  ProjectLDA 

%  Diogo  Pereira  10/22/02 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


clear 
close  all 
clc 

%  W  is  an  array  containing  the  eigenvectors  of  the  covariance  matrix  in  each  column 

%  Load  the  training  images  into  A  matrix 

A=[]; 

T=[]; 

Person=[l  2  3  4  5  6  8  9  1 1  12  13  14  15  16]; 

N_class=length(Person); 

N_object=10; 
for  i=l:N_class 
for  j=l :  10 

im_numl=  num2str(Person(i)); 
im_num2=  num2str(j); 

imgname  =  strcat(im_numl,'-',im_num2,'-l-a'); 
img=  imread(img_name,'bmp'); 

[diml,dhTt2]  =  size(img); 
x=reshape(img,diml  *dim2, 1); 

A=[A  x];  %  A  contains  one  image  in  each  row 
T=[T  Person(i)]; 
end 
end 

clear  x; 

A=double(A); 

N_pictures=N_class*N_object;  %N_pictures  is  the  total  number  of  pictures  on  the  A  ar¬ 
ray 

%show  up  to  100  image  contained  in  the  array  A 
figure(  1 ) 

for  i=l:min([N_pictures  100]); 
subplot(10,10,i) 
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trial_face=reshape(  A( :  ,i),dim  1  ,dim2 ) ; 
imagesc(trial_face) 
colormap(gray); 
axis  off; 
end 

H=diml; 

%  Compute  the  pea  of  the  data  matrix 
%  Amg=mean(mean(A)); 
temp=(N_class*N_object)  -  N_class; 

[W,m,Amean,Ad,EVA]=pca(A,temp);  %  This  reduces  the  dimension  to  N-c 

%  Compute  the  matrix  P  obtained  with  the  PCA 

P=W'*(Ad); 

%  Compute  the  Wopt  of  the  Fisher's  Linear  Discriminant 
C  l=ones(  1  ,N_object); 

C=[]; 

for  temp=l  :N_class 
C=[C  temp*Cl]; 

%C=[C1  2*0  3*0  4*0  5*0  6*0  ...]; 
end 

[Wopt,D]=fld(P,C); 

%  Compute  the  projection  matrix  P  obtained  with  the  FLD(Fisher  Linear  Discriminant) 
P=Wopt'*P; 

%  Computing  the  centroid  of  each  class=the  average  of  the  coefficients 
C2=[]; 

for  i=l  :N_class 

U 1  =(i- 1 )  *N_obj  ect+ 1 ; 

U2=U  1  +N_obj  ect- 1 ; 

C2(:,i)=mean(P(:,(U  1  :U2)),2); 
end 

W2=W*Wopt; 

save  faces2  W2  C2  %Amg 

%  Load  the  test  images  into  B  matrix 

B=[]; 

T=[]; 

Person=[l  2  3  4  5  6  8  9  1 1  12  13  14  15  16]; 

N_class=length(Person); 

N_object=10; 
for  i=l  :N_class 
for  j=l :  10 

im_numl=  num2str(Person(i)); 
im_num2=  num2str(j); 

img_name  =  strcat(im_numl,'-',im_num2,'-5-a'); 
img=  imread(img_name,’bmp'); 

[diml,dim2]  =  size(img); 
x=reshape(img,diml  *dim2, 1); 

B=[B  x];  %  B  contains  one  image  in  each  row 
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T=[T  Person(i)]; 
end 
end 

clear  x; 

B=double(B); 

LB=size(B,2); 

B=B -Amean*  ones(  1  ,LB); 

PB=Wopt' *W'*B; 

%  Computing  the  classification 
dmin=100000000; 
kmax=size(PB,  1); 
for  k=l:kmax 

[D]=classify(PB(  1  :k,:),C2(  1  :k,:),dmin); 
error=T-D; 

n_error(k)=sum(  (error~=0)  ,2); 
end 
figure 

plot(n_error,’-b’) 

xlabel('Dimension  of  the  transformation'); 

ylabel(’Number  of  Errors'); 

title('Classification  errors  using  140  test  faces'); 

grid; 

pause(l) 

end; 

%  %  Computing  the  error  in  the  raw  image 
%  %  Load  the  test  images  into  B  matrix 
%  B=[]; 

%  load  faces  W  H  Amean  C 
% 

%  for  i=l:N_class 
%  for  j =(N_obj  ect  +  1 ) :  1 0 
%  im_numl=  num2str(i); 

%  im_num2=  num2str(j); 

%  img_name  =  strcat(im_numl,'-',im_num2,’-l'); 

%  img=  readpgm8(img_name); 

%  %img=locatehead(img);  %locates  the  head 

%  img=double(img); 

%  [m,n]=size(W); 

%  WI=m/H; 

%  [T 1  ,T2,T3]=facemapini(img,Amean,W,H); 

%  F=facemapA(W,H,  Amean,  img,Tl,T2,T3); 

%  %F=facemap(W,H, Amean, img);  %  returns  F  the  facemap 

%  [M,I]=min(F); 

%  [N,yini]=min(M); 

%  xini=I(yini); 

%  face=img(xini:(xini+H-l),yini:(yini+WI-l));  %  obtain  the  face  on  the  image 
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%  [diml,dim2]  =  size(face); 

%  x=reshape(face,diml  *dim2, 1); 

%  B=[B  x];  %  B  contains  one  image  in  each  row 
%  end 
%  end 
%  clear  x; 

%  B=double(B); 

%  %Nonnalizing  the  image 
%  for  i=l:size(B,2) 

%  B(:,i)=B(:,i)+Amg-mean(B(:,l)); 

%  B(:,i)=B(:,i)/max(B(:,i))*256; 

%  end 

% 

%  temp=N_class*(  1 0-N_obj ect); 

%  B=B-Amean*ones(l,temp); 

%  PB=Wopt'*W'*B; 

% 

%  %  Computing  the  classification 
%  dmin=  100000000; 

%  kmax=size(PB,l); 

%  T=[]; 

%  for  temp=l:N_class 
%  for  temp2=  1 :( 1 0-N_obj  ect) 

%  T=[T  temp]; 

%  end 
%  end 
% 

%  for  k=l:kmax 

%  [D]=classify(PB(  1  :k,:),C2(  1  :k,:),dmin); 

%  error=T-D; 

%  n_error(k)=sum((error~=0),2); 

%  end 
%  figure 

%  plot(n_error,'-b') 

%  xlabel('Dimension  of  the  transformation'); 

%  ylabel('Number  of  Errors'); 

%  title('Classification  errors  using  60  test  faces'); 

%  grid; 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


%  LocateClassifyA  is  the  advanced  version  of  LocateClassify 
%  this  routine  locates  a  face  and  classify  it  according  with  the  PCA  method 
%  but  it  uses  FFT 

%  input:  keyboard  ->name  of  the  file  containing  the  image 

%  file  faces->W  eigenvectors,  H  high  of  the  image  on  the  database,  C  centroids  of 
the  classes 
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%  Amean  mean  of  the  data  used  to  obtain  the  eigenvectors 

%  file  faces2->  W2  projection  directions,  C2  centroids  of  the  classes 

%  output:  number  of  the  class  the  object  was  classified  to 

%  Obs.:  you  need  to  run  the  programs  ProjectPCA  and  ProjectLDA  before  running  this 
%  rotine 

%  to  obtain  the  files  faces  and  faces2 
%  respectively 
%  Diogo  Pereira 
%  last  update:  10/09/2002 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


close  all 

clear 

clc 

%  obtaining  the  name  of  the  image 
name=input('Type  the  name  of  the  file  ->',’s'); 

T=3; 

B=readpgm8(name);  %  read  pgm  files 
[HI,WI]=size(B);%  obtaining  the  size  of  the  image 

%  displays  the  image  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

figure('Units','pixels','Position',[100  100  WI  HI]); 
imagesc(B); 

set(gca, 'Units', 'pixels', 'Position', [0  0  WI  HI]); 

colonnap(gray) 

pause(.8) 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


disp( ['Locating  the  face,  wait...’]) 

%  Locates  and  obtains  the  face%%%%%%%%%%%%%%%%%%%%%%%%%% 

B=locatehead(B);  %locates  the  head 

B=double(B); 

load  faces  W  H  Amean  C 

load  faces2  W2  C2 

[m,n]=size(W); 

WI=m/H; 

[T 1  ,T2,T3]=facemapini(B,  Amean, W,H); 

Q=facemapA(W,H,  Amean, B,T  1  ,T2,T3); 

imagesc(Q); 

colonnap(gray) 

[M,I]=min(Q); 

[N,yini]=min(M); 

xini=I(yini); 

face=B(xini:(xini+H-l),yini:(yini+WI-l));  %  obtain  the  face  on  the  image 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


%  displays  the  face  image%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

figure('Units', 'pixels', 'Position', [400  100  WI  H]); 
imagesc(face) 

set(gca, 'Units', 'pixels', 'Position', [0  0  WI  H]); 
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colonnap(gray) 

%  classify  the  face  imag e%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

face=reshape(face,WI*H,  1); 

face=face-Amean; 

%classify  using  the  PCA  scheme 

PB=W'*face; 

dmin=100000000; 

kmax=size(PB,l); 

k=kmax; 

D=classify(PB(l  :k,:),C(  1  :k,:),dnhn); 

disp(['You  belongs  to  class  ’,num2str(D),'  using  the  PCA  classification  scheme’]) 

%classify  using  the  LDA  scheme 

PC=(W2)'*face; 

dmin=100000000; 

kmax=size(PC,l); 

k=kmax; 

D=classify(PC(  1  :k,  :),C2(  1  :k,:),dmin); 

disp(['You  belongs  to  class  ’,num2str(D),'  using  the  LDA  classification  scheme’]) 
end; 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


%  InfraredTotalExitance 

%  This  program  computes  the  Exitance  of  a  blackbody 
%  plotting  it  within  a  wavelength  range  for  a  given  temperature 


% 


%  Diogo  Pereira 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


clc; 

close  all; 
clear; 

h=6.625e-34; 

k=1.38e-23; 

c=3e8; 

cl=2*pi*h*cA2*  10A20 
c2=c*h/k*10A6 


%  cl=3.741e4; 

%  c2=1.438e4; 

%  The  temperatures  to  be  plotted  (insert  the  desired  temperatures  on  the  array  T 
T=[l  150]; 

LT=length(T); 

LI  =0.5;  %  inferior  limit  for  the  graph  in  micrometers 
L2=30;  %  superior  limit  for  the  graph  in  micrometers 
L=logspace(log  1 0(L  1  ),log  1 0(L2),2000); 

LL=length(L); 
for  p=l:LT 
for  q=l:LL 
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W(q,p)=cl/(  L(q)A5  *(exp(  c2/(L(q)*T(p)))  -  1)); 
end 
end 

%  Generates  a  Matlab  plot 
figure 

plot(L(:),W(:,  1  ),'b-'); 
grid; 

xlabel('Wavelength  [x  10A-6  m]’); 

ylabel(’M  -  Exitance  [Watts/(cmA2  x  10A-6  m)]’); 

title(’Exitance  for  a  blackbody’); 

%  creates  the  legend 
legend('l  150K’); 

%  Save  the  data  on  a  spreadsheet  that  can  be  read  by  excel 
filename- ExitanceBlackBody  1 1 50K’; 

M=[L’W(:,1)]; 

wk  1  write(filename,M,2,2) ; 


function  image  =  readpgm8(filename) 

%READPGM8  Read  a  raw  pgm  file  as  a  matrix 

% 

%  IMAGE  =  READPGM(FILENAME)  reads  the  binary  PGM  image  data  from 
%  the  file  named  FILENAME  and  returns  the  image  as  a  2-dimensional 
%  array  of  integers  IMAGE.  Assumes  the  file  is  a  raw  PGM  file 

%  containing  8-bit  unsigned  character  data  to  represent  pixel  values. 

% 

%  Matthew  Dailey,  1997 
%  Modified  by  Diogo  Pereira,  2002 
%  Open  the  file 
fid  =  fopen(filename,'r’); 

%  Parse  and  check  the  header  infonnation.  No  #  comments  allowed. 

A  =  fgets(fid); 
if  strcmp(A(l:2),'P5')  ~=  1 
error('File  is  not  a  raw  PGM'); 
end; 

A  =  fgets(fid); 
sizes  =  sscanf(A,’%d'); 
w  =  sizes(l); 
h  =  sizes(2); 

A  =  fgets(fid); 
max  =  sscanf(A,'%d’); 
tlength  =  w*h; 
if  max  ~=  255 

error('Cannot  handle  anything  but  8-bit  graymaps'); 
end; 

%  Read  the  raw  data 
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[v, count]  =  fread(fid,inf,'uint8=>uint8');  %  this  makes  v  a  uint8  instead  of  double 
%if  count  ~=  tlength 

%  error('File  size  does  not  agree  with  specified  dimensions.'); 

%end; 

v=v(l:tlength,l); 

%  Pack  the  column  vector  v  into  the  image  matrix 

image  =  reshape(v,w,h)'; 

fclose(fid); 

return 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


%  CropFaceAuto  -This  file  loads  all  the  images  stored  in  pgm  format  and  save 
%  the  croped  image  in  bmp  format 


% 


%  Diogo  Pereira 
%  12/14/2002 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


Person=[l  2345689  11  12  13  14  15  16];%  Person  numbers  contained  on  the  files 

N_pictures=10;  %  number  of  pictures  for  each  person  number 

Section=l; 

load  faces  W  H  Amean  C  %Amg 
for  i=l  dength(Person) 
for  j= 1  :N_pictures 

im  num  I  =  num2str(Person(i));  im_num2=  num2str(j);im_num3=num2str( Section); 
imgname  =  strcat(im_numl,'-',im_num2,'-',im_num3); 

A=readpgm8(img_name);  %  read  pgm  files 
%Head=locatehead(A);  %locates  the  head 
B=double(A); 

[m,n]=size(W); 

WI=m/H; 

[Tl,T2,T3]=facemapini(B,Amean,W,H); 

Q=facemap  A(W,H,  Amean,  B,T  1  ,T2,T3); 

[M,I]=min(Q); 

[N,yini]=min(M); 

xini=I(yini); 

face=A(xini:(xini+H-l),yini:(yini+WI-l));  %  obtain  the  face  on  the  image 
namefinal=[img_name  '-b.bmp']; 
imwrite(face,namefinal,'bmp');  %  save  the  files  as  bmp 
end 
end 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

%  Cross  Validate 

%  This  programs  test  the  2  schemes  using  cross  validation 
%  without  keeping  the  number  of  elements  per  class  constant 
%  Diogo  Pereira 
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%  Load  the  images %%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

A=[];  %  A  contains  images  in  columns 

T=[];  %  T  containt  the  class  number  of  each  picture 

Crop='a'; 

N_pictures=10;  %  number  of  pictures  for  each  person  number 


Person=[l  2345689  11  12  13  14  15  16];%  Person  numbers  contained  on  the  files 
Section=l; 

[A,T]=LoadImage(A,T, Person,  N_pictures, Section, Crop); 

Section=5; 

[A,T]=LoadImage(A,T,Person,  N_pictures, Section, Crop); 

Section=6; 

[A, T]=LoadImage(A,T, Person,  N_pictures, Section, Crop); 

%  Obtain  the  samples  to  be  used  on  the  testing  and  remove  it  from  the  training  images 
N_trials=100; 

N_samples=100; 

APermanent=A; 

TPermanent=T ; 
for  i=l  :N_trials 

B=[];  %  B  contains  the  samples 

TB=[];  %  TB  contains  the  class  number  of  each  sample 

A=APermanent; 

T=TPermanent; 

N_images=length(T); 
for  j= 1  :N_samples 

Sample=round((N_images-  l)*rand)+ 1 ; 

B=[B  A(:, Sample)]; 

A(:,Sample)=[]; 

TB=[TB  T(Sample)]; 

T(Sample)=[]; 

N_images=N_images- 1 ; 
end 


% PCA %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


%  Compute  the  pea  of  the  training  images  A 
[ W,m,  Amean,  Ad]=pca(  A) ; 

%  Compute  the  projection  matrix  P 
P=W'*(Ad); 

%  Compute  the  centroid  of  each  class  contained  in  Person 
C=meanclass(P,T,  Person); 

%  Using  the  testing  images  B 
Ad=B-Amean*ones(  1  ,size(B,2)); 

%  Compute  the  projection  matrix  P 
P=W'*(Ad); 

%  Classify  the  testing  images 
dmin=100000000; 
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[D]=classify(P,C,dmin); 

NZD=find(D~=0); 

D(NZD)=Person(D(NZD)); 

error=TB-D; 

n_errorPCA(i)=sum((error~=0),2); 

%0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o%%% 

%  LDA %%%%%%%%%%%%%%%% 

%  Compute  the  pea  of  the  training  images 
temp=(N_images- 1 )  *  length(Person); 

[W,m,Amean,Ad,EVA]=pca(A,temp);  %  This  reduces  the  dimension  to  N-c  =  36  - 
6=30  the  best  result  was  with  29 
%  Compute  the  projection  matrix  P 
P=W'*(Ad); 

%  Compute  the  projection  matrix  P  obtained  with  the  FLD(Fisher  Linear  Discriminant) 
[Wopt,D]=fld(P,T); 

P=Wopt'*P; 

%  Computing  the  centroid  of  each  class 
C=meanclass(P,T,Person); 

%  Using  the  testing  data 
Ad=B-Amean*ones(  1  ,size(B,2)); 

P=W'*(Ad); 

P=Wopt'*P; 

%  Classify  the  testing  images 
dmin=100000000; 

[D]=classify(P,C,dmin); 

NZD=find(D~=0); 

D(NZD)=Person(D(NZD)); 

error=TB-D; 

n_errorLDA(i)=sum((error~=0 )  ,2) ; 

%%%%%%%%%%%%%%%%%%%%%%%%% 

end 

figure 

L 1  =min(n_errorPCA); 

L2=max(n_errorPCA); 
temp  1 =linspace(0,L2,L2+ 1 ); 

N  l=hist(n_errorPCA,temp  1); 
plot(temp  1  ,N  1  ,’bd-'); 
title('Histogram  of  the  number  of  errors'); 
xlabel(’Number  of  errors'); 
ylabel('Number  of  occurence’); 
hold  on 

L 1  =min(n_errorLD  A); 

L2=max(n_errorLD  A) ; 
temp2=linspace(0,L2,L2+l); 

N  2=hist(n_errorLD  A,  temp2 ) ; 
plot(temp2,N2,’rd-'); 
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legend('PCA  scheme',  'LDA  scheme’); 


%  Function  Facemap  computes  the  square  of  the  difference 
%  between  a  subimage  and  its  projection  on  the  facespace; 

%  the  subimage  has  the  size  of  the  eigenface 

%  input:  A-  array  (m  x  n)  containing  the  eigenfaces  in  columns 

%  m  dimension  of  each  eigenface,  n  number  of  eigenfaces 

%  diml  -  dimension  1  of  each  eigenface,  dim2=m/diml 

%  Amean-  vector  containing  the  mean  of  the  data  used  to  create  the  eigenfaces 

%  B-  array  (p  x  q)  containing  the  image  to  be  analyzed 

%  output  C-  array  (r  x  t)  containing  the  difference  between  the  subimage  and  its 

%  projection  squared 

% 

%  Diogo  Pereira 
%  09/16/02 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


function  [C]=facemap(A, diml , Amean, B) 

[m,n]=size(A); 
dim2=m/diml ; 
tp,q]=size(B); 

%  computing  r  the  number  of  lines  of  our  facemap 
r=p-diml+l; 

%  computing  t  the  number  of  columns  of  our  facemap 
t=q-dim2+l; 

%  creating  an  empty  face  map 
C=zeros(r,t); 

%scanning  the  image 
for  L=l:r 
for  U=l:t 

%  obtaining  the  subimage  PHI 

PHI=B(L:(L+diml-l),U:(U+dim2-l));  %PHI  has  dimension  diml  x  dim2 
%  obtaining  the  subimage  projected  on  the  facespace  PHIT 

PHI=reshape(PHI,  diml*dim2,l)-Amean;  %PHI  is  reshaped  to  (diml  x  dim2)  x 

%1 


Temp=A'*PHI;  %  Temp  contain  the  projection  of  PHI  in  each  eigen¬ 

face 

PHIT=A*Temp;  %  PHIT  contain  the  reconstructed  PHI  using  the  ei¬ 

genface 

%  obtaining  the  difference  between  PHI  and  PHIT 

C(L,U)=sum((PHI-PHIT).A2); 

end 

end 

return 
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